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HYPERSONIC DRONE VEHICLE DESIGN: A MULTIDISCIPLINARY EXPERIENCE 


UCLA’S Advanced Aeronautic Design group focussed their efforts on 
design problems of an unmann ed hypersonic vehicle. It is felt that a scaled 
hypersonic drone is necessary to bridge the gap between present theory on 
hypersonics and the future reality of the National Aerospace Plane (NASP) for 
two reasons: (1) to fulfill a need for experimental data in the hypersonic re- 
gime, and (2) to provide a testbed for the scramjet engine which is to be the 
primary mode of propulsion for the NASP. 

The group concentrated on three areas of great concern to NASP design: 
propulsion, thermal management, and flight systems. Problem solving in 
these areas was directed towards design of the drone with the idea that the same 
design techniques could be applied to the NASP. 

A seventy degree swept double delta wing configuration, developed in the 
70’ s at NASA Langley, was chosen as the aerodynamic and geometric model for 
the drone. This vehicle would be air-launched from a B-l at Mach 0.8 and 
48,000 feet, rocket boosted by two internal engines to Mach 10 and 100,000 
feet, and allowed to cruise under power of the scramjet engine until burnout. It 
would then return to base for an unpowered landing. Preliminary energy calcu- 
lations based upon the flight requirements give the drone a gross launch weight 
of 134,000 lb. and an overall length of 85 feet. 
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1. Introduction 


1 . Introduction 


The efforts of the Advanced Aeronautics Propulsion Design 
group have been directed towards developing tools to analyze the effects 
of aircraft geometry and freestream conditions on scramjet performance. 
Specifically, the group has been concerned with the shock structure, 
boundary layer growth and inlet geometry of a Mach 10 drone research 
vehicle, and the impact of these parameters on the efficiency and thrust of 
the scramjets. 

The tools used to analyze the flowfield have primarily 
consisted of pre-existing computer software. Six programs (STUB, INLET, 
CPIPE, SEAGULL, SCRAM, and EDDYBL) were used to evaluate different 
components of the drone's underbody and scramjet engines. SCRAM is a 
one-dimensional program for the estimation of scramjet performance of a 
hypersonic vehicle. SCRAM begins with the geometry of the vehicle and 
calculates gross estimates of thrust, efficiency, and other engine 
parameters. It is used by our group as a check for the more detailed 
analysis. SEAGULL calculates inviscid conical shock structures and 

freestream conditions for axisymmetric geometries. SEAGULL is used to 
calculate the shock structures over the drone forebody. EDDYBL is a 
two-dimensional/axisymmetric, compressible, turbulent boundary layer 
program, and is used in conjunction with SEAGULL to simulate the viscous 
flow conditions along the forebody. Using the freestream conditions from 
SEAGULL,^ EDDYBL is used to calculate the boundary layer growth along the 
vehicle. INLET uses an inviscid analysis to design hypersonic engine inlets. 
For given*, inlet conditions and desired combustor entrance conditions, 
INLET calculates the cowl geometry required. Data from SEAGULL and 
EDDYBL will be used to generate inlet conditions for use with INLET. 


2. Hypersonic Cruise Vehicle(HCV) 
Configuration 


Finally, CPIPE is a one-dimensional code for the analysis of scramjet 
combustor performance. CPIPE will use the output conditions from INLET 
to calculate the combustor efficiencies. An overall momentum balance on 
the engine should give the actual thrust of the scramjet engine 
configuration. The use of these codes allows a detailed analysis of 
scramjet engine performance from nose to tail of the vehicle. 

An initial application of these codes was performed on a 
five-degree half-angle cone with two five-degree compression ramps 
leading into the inlet. Boundary layer and flow properties were determined 
for the forebody, and an inlet configuration was generated. 

Analysis of this initial design revealed that greater forebody 
compression was required to produce an acceptable inlet design. 
Therefore, the design was altered by placing one additional five and two 
additional three-degree ramps between the nose of the vehicle and the two 
five-degree ramps leading into the inlet. The location of the three-degree 
ramps were varied to maximize static pressure gains. The final design 
provided sufficient thrust to maintain level flight at design condition. 


2. Hypersonic Cruise Vehicle (HCV) Configuration 
It was determined early in Winter Quarter 1988 that the 
Advanced Aeronautics Design Group (AAD) needed a unifying concept for the 
design of the hypersonic drone. This was accomplished by choosing a 
configuration that had been researched extensively at NASA Langley 
Research Center in the late 1970's. The advantage of a single tested 
configuration for the entire class is clear: each group would base their 

designs on the same vehicle, with the same restrictions applying to all 
three groups, thereby forcing greater interaction among the groups. 

The configuration chosen by the AAD students was a 70 degree 
swept delta wing design conceived by NASA Langley engineers. It was a 
hypersonic research airplane concept, and has been dubbed the Hypersonic 
Cruise Vehicle (HCV) by the AAD students. Figure 1 shows the HCV as it 
was originally conceived. The HCV was designed to be air-launched from a 
B-52 at Mach 0.8 and rocket-accelerated up to a cruise speed of Mach 7 by 
four liquid hydrogen rocket motors mounted in the base. At this point, the 
scramjets would be ignited and the vehicle would cruise until burn-out, at 
which point it would return to base for an unpowered landing. A detailed 
description of the HCV’s configuration is given in Figure 2. In developing 
the HCV design, a wind tunnel model was built, and performance 

characteristics of the model were collected and tabulated in Reference 1. 

The geometric characteristics of the model are tabulated in Table 1. The 
scaling for the model was 0.021 of the size conceived for the actual HCV. 

The first step in adapting the HCV for the purpose of a 

hypersonic drone was to size it according to the needs of the class. This 
included calculating the amount of fuel required to go from the air-launch 
speed to cruise speed. From that, size and weight necessary for the fuel 
requirements were determined . In order to do this, it became necessary to 
determine a means of air-launch. The B-52 has size and weight 

constraints consisting of a 24-ft. wing span limit, a 9-foot fuselage 
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height limit, and a gross weight restriction of 70,000 pounds. This limits 
the size of the drone considerably and, in order to reach Mach 10, it was 
found that the weight limit was exceeded. Alternatives were then 
examined. The 747 and B-1 aircraft were determined to be viable options 
for the air-launch. 

Therefore, the HCV is a basis upon which each group can build 
their design. It gives the Flight Systems group the geometry and spacing 
of the vehicle so that they can determine which sensors and how many are 
required in certain places. This would result in a more detailed and more 
exact design. The Thermal Management group can determine critical areas 
on the HCV that would require cooling. They can change the design within 
certain reasonable limits to satisfy their cooling system design. A 

definite configuration for the the drone gives the Propulsion group a 
starting point on the design of the underbelly for the vehicle; that is, a 
first iteration could be performed from which they refined their design to 
achieve an optimum configuration. 
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3. Program Applications 


^.Program Ap plic ations 
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The group explored pre-existing computer programs, as well as 
devising others as needed, as a method of analyzing the merits of various 
designs because of the time savings afforded and the ability to calculate 
complex flowfields that are inherent to hypersonic design. Figure 3 shows 
the areas of application of each of the codes in analyzing the drone. A 
brief description of each program follows, while a more detailed 
description fo the workings of each code may be found in Appendix A. 

SCRAM 

The program SCRAM is a one-dimensional airframe-integrated 
scramjet simulation program. It is currently being used by NASA on the 
National Aerospace Plane (NASP) project as a first iteration on scramjet 

it j.g available to the AAD (Advanced 
Aeronautics Design) students on IBM PC computers and has already been run 
successfully using the HCV geometry as inputs (see Appendix A1 for a 
description of the program and a sample case). 

.The inputs to SCRAM include freestream conditions, scramjet 
geometry, and other parameters necessary for calculating the flowfield 
properties such as kinetic energy efficiency, which approximates 
momentum losses due to shocks, and temperatures along the surfaces. 
SCRAM calculates the flowfield characteristics for five stages of the 
scramjet: freestream, forebody, inlet, combustor, and nozzle. It outputs 

the flow conditions at the end of each stage, as well as boundary layer 
thicknesses. It then calculates the engine cycle performance parameters 
and outputs them in an overall summary. 

SCRAM was used in two ways for the Hypersonic Drone concept. 
First, it was used to help determine the sizing of the engine. By developing 
a scaled drone instead of a full-sized, piloted hypersonic vehicle, there 
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was a concern that the scramjet engine will not be large enough to 
produce sufficient thrust. SCRAM can be run at the smallest size limit of 
the drone to determine feasibility of the design. Since SCRAM generates 
values that are optimistic for the engine performance, it offers a means of 
evaluating relative worth of designs. That is, if SCRAM's output shows 
that a design is not feasible, then the design will be certain not to work, 
and further study into that configuration would probably not be 
worthwhile. 


£IUB 


The flow behind the detached bow shock structure generated by a 
hypersonic blunt body contains regions of both subsonic and supersonic 
flow. Since SEAGULL is restricted to the case of flow that is purely 
supersonic after a shock, there arises a need for an analysis of the shock 
structure and flow properties at the nose of the hypersonic drone studied 
by the Advanced Aeronautics Design class. The shock structure and flow 
properties then become inputs for the SEAGULL and EDDYBL programs. 

STUB is a written program that utilizes the time-dependent 
technique to generate the detached bow shock structure off a hypersonic 
blunt body. Details of the method used by the program can be found in 
References 2, 3, and 4, and the governing equations in Appendix A2. The 
program is also listed in Appendix A2. 

Briefly, the time-dependent technique employs a numerical 
calculation of flow properties through time steps, using the unsteady 
Navier-Stokes equations as the governing equations. Although the 
introduction of time appears to be an added complication for a 
steady-state flow, it becomes a necessary simplification to the problem 
of hypersonic blunt bodies. There currently exists no uniformly valid, 
steady state technique that can handle the mixed supersonic and subsonic 
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region behind the bow shock (Ref. 2). 

Although STUB was written to analyze two-dimensional blunt 
bodies, an extension to axisymmetric or three-dimensional asymmetric 
blunt bodies could be made in the same manner as References 4 and 5. The 
two-dimensional nature of this analysis, however, does mean that 
relieving effects and other three-dimensional, axisymmetric effects for a 
spherically-tipped nose (as should be modelled for the hypersonic drone) 
are neglected. This then became a severe limitation of the program. 


SEAGULL 

SEAGULL is a code designed for the analysis of a two-dimensional 
or axisymmetric supersonic inviscid flow of an ideal gas. For the purposes 
of the Propu lsion Design Group, SEAGULL was used to 

model the external compression along the forebody of a scramjet engine. 
SEAGULL provides only an inviscid solution for the forebody flowfield and 
shock structure, and thus the boundary layer code EDDYBL was used with 
SEAGULL to iteratively solve for a viscous flowfield solution. 

The resulting flowfield properties calculated using SEAGULL (and 
EDDYBL) were used as inputs for the program INLET. One limitation of 
SEAGULL is that it models both an upper and lower wall; thus, the flow 
properties that SEAGULL outputs are for internal compression. This causes 
the flow on the body downstream of the shock to slow, whereas for purely 
external flow, it would remain constant. The slowing effect is minimized, 
however, by placing the walls far apart. 
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EDDYBL 


EDDYBL is a fully compressible, two-dimensional or axisymmetric 
computer program for the calculation of boundary layer properties. The 
program accounts for mass flux at the body surface, heat flux, pressure 
and temperature gradients, and both transverse and longitudinal body 
curvature. The program does not, however, account for the effects of 
shock-boundary layer interaction or separated flow. EDDYBL is currently a 
fully operational program and has been run for supersonic cases. For a 
more complete description of the program EDDYBL and sample output 
parameters, see Appendix A4. 

For the purposes of the Advanced Aeronautics Propulsion Design 
group, EDDYBL was used in conjunction with the program SEAGULL, which 
performed an inviscid analysis on the forebody of the scramjet. Provided 
with upstream conditions and geometry from SEAGULL, EDDYBL output the 
boundary layer characteristics of the flow field. The boundary layer 
thickness was used to redefine the geometry, thus creating a new inviscid 
input geometry for SEAGULL. This will alter the shock structure 
calculated for the forebody, and thus will alter the upstream conditions 
used for EDDYBL. This iterative procedure was used until the flowfield 
properties converged on a viscous solution. 

EDDYBL was helpful in determining the overall performance of a 
given scramjet configuration by calculating the boundary layer growth at 
the inlet entrance (or cowl lip), and thus determining the effective mass 
flow entering the engine ( Fig. 4). Overall, using EDDYBL in conjunction 
with the inviscid codes SEAGULL and INLET will provide a more realistic 
evaluation, of the performance of a given scramjet configuration. 



INLET 


INLET was the code used in designing the internal cowl and 
centerbody geometries of the scramjet. It is a two dimensional inviscid 
code that uses the method of characteristics to determine the internal 
geometry of the inlet (see Appendix A5 for a description of the program 
and sample inputs and outputs). Given the conditions at the entrance to the 
inlet and the desired combustor entrance conditions, INLET calculates the 
cowl and centerbody geometries required to minimize total pressure loss 
across the inlet. Although the distinction between the forebody and inlet 
can be nebulous, the propulsion group decided to use SEAGULL mainly for 
external compression and INLET for internal- compression. Thus, the main 
benefit of using INLET is that it determines the cowl and centerbody 
geometries for the scramjet engine. However, it was necessary to meet 
the desired combustor entrance requirements, as well as the upstream 
entrance conditions (acquired from SEAGULL). Therefore, it was necessary 
to iteratively use INLET and CPIPE (combustor code) to determine whether 
INLET could provide reasonable cowl and centerbody geometries for the 
combustor requirements. 


CPIPE 


CPIPE is a one dimensional real gas analysis of the combustion of 
hydrogen in air in a supersonic channel. Given the initial upstream 
conditions and geometry definition, CPIPE outputs the downstream one 
dimensional flow properties (see Appendix A6 for a description of the 
program and required inputs). 

CPIPE was used to model the combustor region of the scramjet 
engine. CPIPE is a relatively simple analysis of the combustor since it 
uses only one dimensional conservation equations to calculate the 
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downstream conditions. However, more complicated analyses of the 
supersonic combustor are rare and tend to concentrate on the flow 
structure in specific areas of the combustor (i.e. directly behind a 
flameholding structure or fuel injection area). In addition, CPIPE has 
several useful features such as (1) possibilities for several fuel injection 
points at a variety of angles, (2) estimates of ignition and flameholding 
• likelihood, and (3) estimates of heat flux, given temperature distribution. 
Thus, CPIPE served the purpose of scramjet performance analysis without 
being unnecessarily complex. 

Since CPIPE is one dimensional it was necessary to average the 
flowfield values obtained from INLET in order to input the appropriate 
values in CPIPE. 

References 10, 11, and 12, which contain temperature and 
pressure distributions for various combustor designs that were tested in 
Langley's supersonic combustor test facility, could be used to validate 
CPIPE. This valdation was not completed due to time constraints! however, 
CPIPE has been used satisfactorily at NASA Langley to simulate scramjet 
combustor designs. 
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SCRAM was used as a first iteration on the scramjet design for 
the drone. It was primarily used to determine whether an engine small 
enough to fit onto the HCV would produce enough thrust to overcome the 
total aircraft drag. One of the primary weaknesses of SCRAM, however, is 
it's dependance on wise choices of abstract efficiency parameters. Such 
parameters define the overall kinetic energy efficiency of the design, the 
pressure loss coefficients across the various shocks and through the 
combustors, and have a large impact upon the overall thrust delivered. 

When SCRAM was run with the HCV's geometry and flight 
condition as input, it was found that the scramjet could supply the 
required thrust. With a nozzle efficiency of 1.0, SCRAM predicted a total 
thrust force of 43,000 pounds. This value was as low as 30,000 pounds 
with a nozzle efficiency of 0.96. Using information given from wind tunnel 
tests on the HCV configuration, [13], the total drag at Mach 10 and 100,000 
feet was found to be 26,600 pounds, for an excess thrust of 3,400 pounds 
with 96% nozzle efficiency. 

The vehicle nose was modelled as a two-dimensional, five degree 
ramp, with a rounded, one-inch radius nosetip. Initial guesses for the 
shock structure and flow properties were inputted into STUB, and final 
converged values were computed. These values for the shock structure and 
flow properties were made available too late for use, but would have 
served as inputs for SEAGULL and EDDYBL. As it was, the vehicle was 
modeled as having a sharp nose, compared to the one-inch radius nose of 
the actual vehicle. The next step, then, would have been to proceed through 
one iteration using STUB outputs. 

The group began to analyze the merits of various forebody geometries. 
For the following forebody geometries presented, two opposing constraints 
were considered. Static compression of the flow by the forebody had to be 
attained with maximum efficiency and minimal total pressure loss. At the 



same time, the flow had to be slowed sufficiently prior to the inlet so that 
the cowl would be made as short as possible. 

The first geometry analyzed was a simple five degree cone and also 
acted as verification for SEAGULL. (Appendix A3). The cone did not 
compress or slow the flow sufficiently and required a cowl length of 150 
feet to compress the flow from Mach 9 ahead of the inlet ramps to the 
maximum combustor speed of Mach 4. 

The second geometry consisted of two ramps in addition to the basic 
five degree cone. (Fig. 5). The first ramp is five degrees and is located 
fifteen feet from the nose. A second, three degree ramp is located 
twenty-four feet from the nose for a total turning angle of thirteen 
degrees. This design resulted in a total pressure loss of 22%, and a static 
pressure increase of 8.97 before the inlet ramps. The Mach at the entrance 
to the inlet was 7.51. (Table 2). 

The third geometry moved the ramps of the second design three feet 

m rnmmmmm mrn n - 

further from the nose (Fig 6). This resulted in a decrease of total pressure 
loss at the inlet and an increase in static pressure. (Table 2). These 
changes improved the performance of the forebody as a compression 
surface. The last geometry added another three degree ramp thirty-one 
feet from the nose. This resulted in a total pressure loss of 23%, but a 
static pressure increase of 14.6. (Table 1). This geometry was the one 
chosen for the final design. (Fig 7). 

Boundary layer properties were computed for a variety of geometries 
in the early stages of the design. The simplest case run was for the five 
degree cone which was modeled as having a sharp nose, unlike the rounded 
nose of the drone. This approximation introduces error in the downstream 
flowfield, but was necessary because of a lack of freestream data in that 
region. The surface rougness of the vehicle was set to .001 in., the default 
value of EDDYBL, which is a very smooth surface. Lacking any data for high 
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altitude atmospheric turbulence, the freesteam turbulence intensity was 
set to 2 %, which is what could be expected from a well designed wind 
tunnel [14]. 

Tha results of the computations are presented in Figures 8 and 9. Note 
first the "blips" that occur for the boundary layer thickness. These are 
points where the program converged on an incorrect thickness, but with 
the next step corrected itself and converged on the correct value. 
Communication with the author of EDDYBL indicates that these points are 
stations that have difficulty converging due to the grid representation of 
the flowfield which result in numerical errors. The automatic addition of 

grid points by the program can cause EDDYBL not to converge -for that 
station, but generally converges to the correct solution with the next 
step. 

Note also the trend of the lower left corners of the curves to have a 
slight lip. This is also caused by convergence errors of the program. 
EDDYBL calculates the stagnation properties from input free stream 
conditions and uses these initial starting conditions to predict the 
boundary layer properties at the next step. If that prediction is incorrect 
and the program does not converge after twenty iterations, the program 
will proceed to the next station. To say that the program does not converge 
does not mean that the calculations diverged, merely that the solution did 
not converge fast enough. After several stations, however, the program 
begins to make to make correct predictions and the solution converges. 
Extrapolation of the curves, shown by dotted lines in Figures 8 and 9 shows 
that the computed boundary layer values are close to zero at the stagnation 
point, as they should be. 

The. shape of the curves in Figure 9 suggests that transition does not 
occur and that the flow is entirely laminar. Literature shows that 
transition on 10 degree cones occurs at an arc-length Reynolds number of 
approximately ten million, while EDDYBL shows that transition has not yet 
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occu r od at Reynolds numbers of 48 million^ 15, pp 39, 59,61]. The data 
presented In the literature, however, is for cones with roughness heights 
of .015 in., and unknown freestream turbulence intensity, both of which 

have an impact on transition to turbulence and are different than those 
used for this case. 

Figure 10 shows the heat flux at the surface of the vehicle as a 
function of axial distance, with the surface maintained at 1000 degrees 
Rankine. ( A misunderstanding occured between the thermal management 
and propulsion groups as to the temperatuure that the wall would be 
maintained at. The thermal management group specified 1000 degrees 

Kelvin, which is higher than the 1000 degrees Rankine used here. The heat 
flux could thus be expected to be somewhat lower than that indicated.) The 
portion of the curve to the left of point A of Figure 10 is incorrect, 
suffering from the same convergence problems near the stagnation point as 
the boundary layer properties of Figures 8 and 9. The heating rates shown 
have been determined to be acceptable by the thermal management group. 
Active cooling would be required to maintain the skin temperature but 
would be within the capabilities of the cooling system. 

Figures 11 through 14 show the effect of boundary layer suction at 
the wall. The mass removal began 31 feet from the nose at a rate of .0005 
s!ugs/ft A 2*s. (Figure 13). Figures 11 and 12 show the effects of the mass 
removal on boundary layer thickness. The impact is seen to be primarily on 
displacement, and is caused because the fluid that is being removed is the 
slowest and has the largest impact on displacement thickness. (Recall that 
displacement thickness is a measure of the mass flow rate along the body 
that is lost due to the slowing of the flow in the boundary layer.) 

The mass removal has a profound effect on heat flux. Figure 14 shows 
that, although the suction has a small effect on boundary layer thickness, 
the heating rate at the wall was increased to a value comparable to the 
heating rates experienced in the first few feet of the forebody. Increased 
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boundary layer suction might therefore lead to unacceptably high heating 
rates, which would put an upper limit on the amount of boundary layer 
control that could be obtained through suction. 

Finally, several hand-calculations were performed to estimate the 
total mass removed. (Appendix B). The results indicate that, for the five 
degree cone, the mass removal rate per unit width over the last five feet 
of the forebody directly before the inlet ramps would be .00425 
slugs/sec*ft. Assuming a depth of 7 feet, (the width of the forebody at the 
inlet ramps), the total mass flow rate would be .0298 slugs/sec. The 
average velocity of the removed mass would be 57.2 ft/sec, thus removing 

400 ft A 3/sec over the 35 square foot area, which would be a large volume 
of air to dispose of. 

Another calculation was made to determine the removal rate that 
would be required to remove the same amount of boundary layer as 
displacement thickness over 7 square feet of the forebody. The required 
mass flow rate would be .208 slugs/sec, with an average velocity of 
400 ft/sec, and a total volume cf 2795 ft A 3/sec. Such a high mass removal 
rate would be unacceptably large. 

At the outset of the Spring quarter, it was still the intention of 
the Propulsion group to use INLET as an axisymmetric code when designing 
and analyzing the scramjet inlet geometry. However, an error in the code 
prevented the group from utilizing INLET'S axisymmetric feature. This 
turned out to benefit the group since the actual nature of each scramjet 
module is more closely represented by a two-dimensional geometry 
because of its location far aft of the nose (see Figure 2) . Once the 
attention of the group was turned to a two-dimensional analysis of the 
scramjet modules, a test case was run using an inlet Mach number of 2.5, 
two ramps of five degrees each, and an exit Mach number of 1.3. The ouput 
of this result was the same as that calculated using the shock relations 
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5. Final Design 



from Reference [ 2]. With the program running correctly, W- first inlet 
geometry case was run. 

The first inlet geometry case was actually a dual case. Since 
problems; were experienced running other computer codes, it was thought 
that the flow could only be slowed to Mach 9.3 at the starting point where 
INLET calculations would begin. Hence, INLET was run for two cases: a 
Mach 9.3 to Mach 3 case and a Mach 9.3 to Mach 4 case. The results of 
these cases are shown in Figures 15 to 22. The output generated was not 
favorable. Inlet plane area to exit plane area ratios were 0.075 for the 
Mach 9.3 to Mach 3 case and 0.173 for the Mach 9.3 to Mach 4 case. This 
means that for an inlet area of 2 feet, the outlet area would be 1.8 inches 
for the Mach 3 exit condition and 4.14 inches for the Mach 4 exit condition. 
These small areas would not allow enough flow through the combustor to 

produce the thrust required to propel the drone at Mach 10 for 5 minutes. 

Not only For these ‘cases a problem, but the total pressure 

losses were enormous due to normal shock waves within the inlet. The 
total pressure recovery for the Mach 9.3 to 4 case was 0.09 and 0.02 for 

the Mach 9.3 to Mach 3 case mainly because of the normal shock wave 
within the inlet. From these results, it was apparent that the free stream 
Mach number needed to be reduced to obtain a reasonable scramjet engine 
design. 

The solution to the problem was to recompute upstream Mach 
numbers. After correcting an operational error with SEAGULL, a free 
stream Mach number of 7.0 could be obtained in order to use as input to 

INLET. The final inlet geometry was generated for a Mach 7 to Mach 4 

compression. The results are presented in the design section of this 
report. 
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The forebody geometry consists of three ramps: one five degree and 
two three degree ramps, as previously mentioned and shown in Figure 7. 
This resulted in 23% loss of total pressure but an increase of 14.23 in 
static pressure. Furthermore, the shocks lie sufficiently close to the body 
so that the inlet capture area is 7.1 times the width of the cowl inlet, and 
there is no spillage for on-design conditions. (Figure 23). 

Difficulties encountered in EDDYBL resulted from the lack of 
shock-boundary layer modeling in the program. (Figures 24 and 25). The 
problem is that the boundary layer properties are discontinuous across the 
shock and it is difficult to determine the correct starting conditions to 
compute properties downstream of the shock. As seen from Figures 25 and 
26, the arc-length Reynolds number increases discontinuously across the 
shock, despite a discontinuous decrease in freestream Mach number. 
However.^aj^ tftffi | second shock, .the Reynolds number decreases with a 
decrease in Mach number, as it should. The result of these variations is 
reflected in the boundary layer properties. The treod in compressible 
boundary layers is for increasing boundary layer thickness with increasing 
Mach number. [2, p. 538]. This agrees with the first discontinuity of Figure 
24, but not with the second. Furthermore, shock-boundary layer interaction 
is not modeled in EDDYBL. Although the computed flow is laminar, 
experiments have shown, very generally speaking, that turbulent boundary 
layers at Mach 7 experience a 10% increase in thickness across a shock. 
[16, p. 7501 -A-1], Note that the computed boundary layer decreased in 
thickness across the shock, showing that EDDYBL is not applicable to 
shock-boundary layer interactions. (Figure 24). 

With the program unable to predict boundary layer properties across 
shocks, quantitative analysis of the design is impossible. Qualitatively, 



however, a number of results may be postulated based upon the results of 
the single five degree cone case. 

Knowing that the boundary layer is .15 ft thick at the first ramp 
location downstream of the nose , boundary layer control will probably be 
required for efficient engine operation. In addition, the boundary layer will 
probably undergo transition prior to entering the inlet, due to traversing 
three shocks and increasing pressure. Forestalling transition, either by 
. boundary layer diverson or suction will therefore be a necessity due to 
both heating and engine performance considerations. 

Control of the boundary layer presents a major problem. Boundary 
layer diversion, as opposed to boundary layer suction, was discarded due 
. to the high cooling expected to be required on the sharp edge of the splitter 
plate. Experience with the five degree cone, however, indicates that large 
cooling rates will be required for boundary layer suctiuon over a very large 
area due to the large removal rates required. This may make splitter plates 
more practical than boundary layer diversion because of it's simplicity. 
Both methods of control would require high cooling rates, as well as a 
large rate and volume of fluid diverted or ingested. The splitter plate 
would accomplish this task without pumps and other active removal 
devices required in boundary layer suction. 

The final configuration for the inlet geometry of the scramjet was 
obtained by running the program INLET for a two-dimensional case by 
specifying a Mach number of 7.0 for the "free-stream" velocity and a Mach 
nubmer if 4.0 for the exit velocity, "free-stream" meaning the flow 
conditions after the initial ramp compressions on the forward underbody of 
the drone but before the two five degree ramps on which INLET constructs 
the flow field. Three internal isentropic compression sections were used 
and the ouptut from INLET is summarized in Figure 27 thru 32. 

Figure 27, Inlet Area vs. Horizontal Distance, shows the inlet 
cross-sectional area measured axially from the cowl lip postion. If the 
graph is scaled by setting the initial area equal to 1.07 square feet, 
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assuming a 1.0 foot depth because it is a two-dimensional case, rnen 0.36 
on the nondimensional scale is 1 .07 feet on the actual drone. This' 1 .07 feet 
dimension is the approximate inlet length shown in Figure 2. This scaling 
will allow the discussion of the results to be related to the actual 

dimensions of the drone. Using the scale 0.36 *= 1.07 feet, it can be seen in 

Figure 27 that the inlet area contracts rapidly from 1.07 feet (0.36) to 

0.456 feet (0.154) in a distance of 3.2 feet (1.08). This sharp decrease in 
area then tapers off and the final inlet area is 0..394 feet (0.133) which 
occurs 20.0 feet (6.75) aft of the initial spike tip. Note that the area is 
approxiately 0.394 square feet (0.133) from a distance of 19.1 (6.44) to 
20.0 (6.75) feet. .This 0.9 feet is needed to ensure uniform properties 
across the exit plane give an area contraction of 2.7 times, meaning that 

the area at the entrance plane is 2.7 times greater than the exit plane area. 

The relative locations of the five degree ramps and the cowl lip 
can be seen in Figure 28. Again, using a scale of 0.36 - 1.07 feet it can be 
seen that the intersections of the two five degree ramps is at 6.87 feet 
(2.32) behind the spike tip of the first ramp and the cowl lip is located 
14.2 feet (4.78) aft of the spike tip of the first ramp. The intersection of 
the two five degree ramps was calculated by INLET so that the shock 
waves generated by the two ramps intersect at the cowl lip. Placing both 
shock waves at the cowl lip minimizes flow spillage and increases engine 
efficiency. Figure 28 also shows that the minimum cross-sectional area of 
the inlet occurs at 18.4 feet (6.2) along the axis of the drone and remains 
constant up to the exit plane at 10.0 feet (6.75). 

Mach number and pressure distributions on the cowl contour are 
shown in Figures 29 and 30, respectively. The magnitude of the Mach 
number can be traced in the output from INLET from the spike tip of the 
first five degree ramp to the exit plane of the inlet. The "free-stream" 
Mach number is 7.0. Behind the shock wace generated from the first five 
degree ramp the Mach number is 6.106. The Mach number after the second 
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es shock wave originating at the intersection of the two five degree- ramps is 
5.401. These results can be calculated using the Prandtl-Meyer oblique 
shock relations found in reference [2]. The shock wave originating from the 
cowl lip further reduces the Mach number to 4.267. This Mach number 
remains constant because INLET assumes that the cowl contour is a 
straight wall for the first 1.73 feet (0.583) for the geometry calculated. 

At 15.88 feet (5.359), the cowl contours change from being a straight wall 
to being an isentropic compression surface and assumes an appropriately 
curved shape. The flow is compressed from a Mach number of 4.267 at 
15.88 feet (5.359) to a Mach number of 4.00 at 17.63 feet (5.950). From 
Figure 30 it can be seen that the pressure increases as Mach number 
decreases. The static pressure ratio increases from 14.9 at the entrance of 
the inlet to 21.2 at the exit plane of the inlet which means that 
approximately one-third of the compression is taking place internally. 

An extremely important feature of all supersonic flow engines is 
the total pressure loss, which dictates the efficiency of the engine. High 
total pressure loss indicates a low engine efficiency while a low total 
pressure loss indicates high engine efficiency. The total pressure ratio 
along the centerbody contour is plotted in Figure 31 and the total pressure 
ratio along the cowl contour is plotted in Figure 32. It can be seen that the 
total pressure ratio along the centerbody contour behind the first shock 
drops to 95% of its initial value behind the second shock. Past the second 
shock, the flow along the centerbody contour encounters no other shock 
waves because the flow is compressed isentropically, which means no 
further total pressure loss occurs. Behind the oblique shock wave at the 
cowl lip the total pressure is only 77% of it’s "free-stream" value. Because 
the cowl lip shock wave is cancelled at the centerbody by INLET, there are 
no internal shock reflections. Since the flow is compressed isentropically 
internally, the total pressure along the cowl surface is constant. Since the 
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to tat. tp: >ure ratio along both the centerbody and cowl contours js 77% at 

the exit plane, the overall total pressure recovery for the scramjet from 
Mach 7 at the spike tip of the first ramp to Mach 4 of the exit plane is 77% 
for the geometry generated. 

The configuration used for the combustor is illustrated in Figure 
33. It is a two. dimensional box combustor with an entrance height of 5 
inches, and a constant width of 11.2 inches ( the combustor length is 10.0 
feet). There are two fuel injection points, at 3.75 ft. and 5.0 ft. 
respectively in the axial direction. Information regarding the injectors is 
given in Table 3. The temperature distribution as a function of axial 
distance is given in Figure 34. There is a steady rise in temperature until 
the location of the first injector, after which it levels off and even begins 
to decrease. This corresponds to a leveling off of the combustor efficiency 
(see Figure 35) at a value near 0.95 after the first injector. This gradual 
rise in colfteti^ttVi^feYfeciency starting at the entrance to the combustor is 
due to the manner in which CPIPE "distributes" the reacted fuel 
equivalence ratio over the combustor. Thus, by the time the flow reaches 
the first injector, a large fraction of the fuel injected at that point has 
started combustion. After the first injector, there is a constant level of 
reaction. This explains the drop in static pressure ( Figure 36) after the 
first fuel injection point. The pressure due to the combustion process is 
constant while the area of the combustor continues to increase. The Mach 
number displays a similar behavior (Figure 37), stabilizing at a constant 
value after the first injector. The exit Mach from the combustor is 
approximately 2.3. 

With the combustor exit properties analyzed by the program 
CPIPE, the flow field inside the nozzle was modeled using the method of 
characteristics. Reference [2] provides the algorithm used to implement 
the computer code. In applying the method of characteristics, the flow 


field is assui -d to be irrotational, inviscid and adiabatic. In addition, the 
characteristic lines are determined to a first order appoximation. 

Once the nozzle exit velocity was determined to be 11,000 
ft/sec, a momentum balance was performed to determine the thrust. The 
control volume used for the momentum balance was defined by the 
equivalent inlet capture area under the nose of the HCV and the vertical 
nozzle exit plane. There was no spillage for the on-design condition; 
hence, it was not included in the momentum balance. 

The thrust available with the nozzle design illustrated in 
Figure 28 was calculated to be 7646 lbs. per unit span, which is 
approximately 8.99% of the thrust available in isentropic expansion. 
Enthalpy losses due to heat flux and viscous effects at the upper nozzle 
wall, as well as plume drag, which was not accounted for, are expected to 
further reduce the efficiency of the nozzle although there should still be 
adequate thrust. 

Note that a certain errors are propagated in locating the mesh 
points with the first order approximation. Reference 17 provides the 
thrust coefficient measured on a similar model tested at Mach 10. The 
thrust coefficient for this design is calculated to be 0.36, which falls 

between 0.1 and 0.92, the range of thrust coefficients presented in 
reference 2. 
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6. Conclusion 


6, Conclusion 


Various forebody geometries were evaluated for their ability to 
efficiently compress free stream flow prior to ingestion by scramjet 
inlets. Also, boundary layer growth on these forebodies was investigated, 
and boundary layer control, in the form of surface suction, was found to be 
somewhat effective. Inlet geometry design was investigated, and it was 
found that inlet speeds greater than Mach 7 would not allow sufficient 
capture area for the compression required. 

A final design for the forebody consisted of a five-degree 
half-angle cone with a five degree compression ramp placed 18 feet aft of 
the nose, and two additional compression ramps placed at 27 and 31 feet. 
The Mach number entering the inlet was 7.2, and an inlet geometry was 
designed using two pre-inlet five degree compression ramps. The cowl 
geometry was generated to yield an isentropic compression to Mach 4 at 
the entrance to the combustor. The total length of the two pre-inlet ramps 
and cowl was twenty feet. 

Although EDDYBL was not able to compute shock-boundary layer 
interaction and was therefore unable to calculate flow properties 
accurately at the shocks generated at the three degree ramps, results from 
the single five degree half-cone were used to obtain trends in boundary 
layer growth on the forebody. Boundary layer suction was found to require 
high suction velocities and large volumes of fluid ingested, accompanied by 
severe increases in surface heating. 

CPIPE was used to compute combustor performance of the 
combustor, and the method of characteristics was used to roughly compute 
the expansion in the nozzle. Finally, a momentum balance was performed 
and the thrust generated by the design was calculated to be 58,000 lb, 
which was slightly higher than the values output by SCRAM, and was 
sufficient to propel the drone at a Mach 10 cruise speed. 
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TABLE I.- GEOMETRIC CHARACTERISTICS OF MODEL 


Uing: 

Area, reference (includes fuselage intercept), m 2 (in 2 ) . . . 0.043 (67.200) 

Area, exposed, m 2 (in 2 ) 0.023 J36.121) 

Area, wetted, m 2 (in 2 ) 0.047 ( 7 2 * 2 J 2 } 

Span, m (in.) O- 21 ? (8 ; 5 * 2 > 

Aspect ratio v ; ; * * _ 1 *^ 

Root chord, at fuselage center line, m (in.) °-353 U3.o9b; 

lip chord, m (in.) \ . 0.085 (3.355 

Taper ratio • • * °.2 1 

Mean aerodynamic chord, m(in.) 0.248 (9.779 } 

Sweepback angles: 

Leading edge, deg ^0 

25- percent- chord line, deg 84 

Trail — — — — — ,? 

Dihedral angle 1 , at airfoil mean line, deg -3.64 

Tnrl'* "M a nfl i Wm ii iVW i m i r i ® 

Airfoil section (See fig. 5(a)) 


Airfoil thickness ratio: 

• Exposed root 

Tip 

Leading-edge radius at - 

Fuselage- line chord, m (in.) 
Tip, m (in.) ....... 

Area of both elevons, nr (in^) 


“ 0.05 

0.06 

5.08 x 10“J (0.020) 

5.08 x 10” M (0.020) 

. . . 0.005 (7.161) 


Forward delta wing: 

Area exposed, outside of fuselage, forward of wing 

leading edge, m 2 (in 2 ) 0.002 (3-394) 

Leading-edge sweep, deg 80 


Tip fin: 

Area, each, m 2 (in 2 ) 

Span, m (in.) 

Aspect ratio 

Root chord, m (in.) 

Tip chord, m(in.) 

Taper ratio 

Mean aerodynamic chord, m (in.) . 


0.0 7 ,.} 


0.004 (5.848) 

0.069 (2.730) 

1.274 

0.086 (3.383) 

0.029 (1.135) 

0.336 

0.062 (2.445) 
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TABLE I.- Concluded 




Sweepback angles: 

Leading edge, top, deg . 
Leading edge, bottom, deg 
Trailing edge, top, deg . 

Toe-in angle, deg 

Airfoil section: 


Leading-edge radius, n (in.) 5.08 * 10"^ (0.020) 

Center vertical tail: 

Area, exposed, m 2 (in 2 ) 0.007 (11.492) 

Span, exposed, d (in.) 0.086 (3^380 ) 

Aspect ratio of exposed area 0.994 

Root chord, at fuselage surface line, n(in.) 0.128 (5.040) 

Tip chord, o (in.) 0.045 (1.760) 

Taper ratio 0.349 

Mean aerodynamic chord of exposed area, a (in.) 0.093 (3.664) 

Sweepback angles: 

Leading edge, deg 55 >0 

Trailing edge, deg * * * * * 24^6 

Airfoil section: 

Thickness ratio at - 

J ip . * ‘ ‘ 0.106 

Root 0.106 

Leading-edge radius, a (in.) 5. 08 * 10" 1 * (0.020) 

Fuselage: 

Length, m (in.) 0.508 (20.000) 

Maximum height, m(in.) 0.071 (2.782) 

Maximum width, m (in.) 0.073 (2.866) 

Fineness ratio of equivalent round body 6.822 

Planform area, m 2 (in 2 ) 0.026 (40.445) 

Wetted area, m 2 (in 2 ) O.O83 (128.460) 

Wetted area, with wing on, m 2 (in 2 ) 0.078 (120.695) 

Wetted area, with both delta wings on, m 2 (in 2 ) . . 0.077 (118.747) 

Base area, m 2 (in 2 ) o!o02 (3!726) 

Complete model, with both delta wings: 

Planform area, m 2 (in 2 ) • , 0.052 (79.960) 

Aspect ratio of planform 0913 


Itr or-, yr A j t'C l - 6 .C 2 / 
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FIGURE 4- Effective Mass Flow entering engine. 










FIGURE 7: Final Forebody Geometry 



FIGURE 8: Boundary Layer Growth vs. Axial Distance on Five Degree Cone 
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FIGURE 9: Boundary Layer Growth vs. Reynolds Number on Five Degree Cone 
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FIGURE 10: Heat Flux on Five Degree Cone 
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FIGURE 1 1 : Boundary Layer Growth vs. Axial Distance on Five Degree Cone 

With Suction 



o 





•FIGURE 12: Boundary Layer Growth vs. Reynolds Number on Five Degree Cone 

With Suction 
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FIGURE 13: Mass Flux on Five Degree Cone 
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FIGURE 17 
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FIGURE 20 
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FIGURE 21 

Mach Number vs. Distance 
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FIGURE 23: Shock Structure on Final Design 
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FIGURE 28 
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FIGURE 29 
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TABLE 3 Information on Injectors 


Injector #1 Injector #2 


Axial 

Distance 

3.75 ft. 

5.0 ft 

Fuel stagnation 

Temperature 

• 

1000 K 

1000 K 

Fuel stagnation 
Pressure 

4758000.0 N/m2 

4758000.0 N/m2 



• 

Mach Number 

1.0 

3.0 

Injector angle 
from horizontal 

0.78 Radians 

0.0 Radians 

Wall angle 
at injector 

0.0 Radians 

1.57 Radians 

t ■ 
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FIGURE 35: Combustor Efficiency vs. Combustor Axial Distance 








FIGURE 36: Combustor Pressure vs. Combustor Axial Distance 
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FIGURE 38 - ISOLATED LIFT ORAO AND MOMENT COCmCICNTS 
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Appendix A1 : SCRAM 



Description 

SCRAM is a one-dimensional airframe-integrated scramjet 
engine simulation program. It was written at the NASA Langley Research 
Center -to give a first-order estimate of scramjet performance. SCRAM 
performs a real gas analysis for a five-stage model. These stages include 
the freestream, the vehicle forebody, the inlet, the combustor, and the 
exhaust nozzle. Once the flow for each stage is analyzed, the overall 
performance is calculated and a summary is given. 

SCRAM is capable of analyzing: (1) oxygen, hydrogen, and air 
injection into the inlet and combustor; (2) axisymmetric or planar 
forebodies; (3) mixed combustion; (4) the boundary layer for a flat plate; 
(5) multistage perpendicular or parallel injection; (6) equilibrium or 
frozen flow referenced by area ratio in the combustor; (7) several internal 
hydrogen fuel mixing models; (8) a swept, symmetric wedge inlet, 
calculating inviscid real gas performance. The program only works for 
started inlets, and flow must be supersonic through the inlet and nozzle 
portions of the code. In addition, SCRAM is not usable if the combustor 
flow becomes subsonic. It is currently configured to handle 
two-dimensional engines only. 

The freestream flowfield is calculated in Subroutine FSC1, 
where, given altitude, Mach number, and other parameters, it computes 
static temperature and pressure from a table based on the 1979 Standard 
Atmospheric Table. The other freestream flow properties are calculated 
from subroutines which find the real gas composition and the dynamic 
state variables. 

The forebody flow properties are calculated in Subroutine 
FLOWF2 : This subroutine accounts for forebody shocks, forebody boundary 
layer heat transfer, and friction losses, given the freestream conditions 
and forebody parameters. It utilizes a specified forebody kinetic energy 
efficiency that approximates engine capture flow momentum losses due to 
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the forebody shock system. Boundary layer properties and heat transfer 
through the boundary layer are calculated iteratively. 

Subroutine INLET1 computes the flow properties for the inlet. 
It models flow through the inlet as supersonic flow undergoing a 
one-dimensional compression process with heat transfer and friction 
losses. Area ratio, length, and kinetic energy efficiency are the required 
inputs. INLET1 also is capable of handling oxygen/hydrogen mass flow 
injection. • 

The flow through the combustor is analyzed in Subroutine 
COMBUST. COMBUST uses perpendicular hydrogen fuel injection for up to 
seven injector stations; one of these stations can also be set for parallel 
injection. It approximates shock losses with the use of the kinetic energy 
efficiencies input by the user. Flow is initially assumed to be in 

equilibrium, but provisions can be made for frozen flow to occur. 
Regenerative cooling may also be taken into account. 

Subroutine NOZZLE1 describes a one-dimensional expansion 
process simulating gas expansion for a scramjet nozzle. It .computes 
boundary layer properties for freestream Mach numbers greater than 6.0, 
but neglects them for freestream Mach numbers less than 6.0. 

Finally, subroutine SCRAP1 computes engine cycle performance 
parameters. It utilizes conservation of momentum equations applied to 
the engine as a whole, using engine geometry and flow properties 
computed by the component routines as inputs. It outputs thrust, specific 
impulse, and the coefficient of thrust as the scramjet performance 
parameters. 
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INTRODUCTION 


Since steady supersonic flow fields are governed by hyperbolic 
differential equations, whereas steady subsonic flow fields are described 
by elliptic differential equations (Reference 2), the mixed region behind 
the detached bow shock structure cannot be handled by any steady-flow 
technique. However, the unsteady Navier-Stokes conservation equations 
become hyperbolic with time, irregardless of the sonic nature of the flow. 
Therefore, the time-dependent technique becomes a simplified means 
through which the steady-state solution is attained. 

E r 0 q r a n i nBe&ctip i ■ 

Details of the time-dependent technique are given in References 
2, 3, and 4, and only a brief description is outlined here. 

The flow field was covered with a mesh and transformed into a 
system of rectangular coordinates. Initial guesses for the shock structure 
and flow properties were inputted and non-dimensionalized. The flow field 
was divided into four regions and treated separately: [1] the shock, treated 
as a discontinuity: [2] the body, examined in the same manner as the shock; 
[3] the inner points between the shock and body, analyzed by the unsteady 
conservation equations: and [4] the upper boundary of the inner points, 
computed through linear extrapolation and restricted to lie beyond the 
sonic line. Time steps were taken in accordance to the CFL stability 
criterion outlined in References 2 and 3. 

First, the shock was analyzed for each time step. At a particular 
time step, an initial guess for the shock velocity was made (the shock will 



"move" through time until it converges onto its "steady-state" location), 
and the new shock location was found. Flow properties at the shock were 
determined by the Rankine-Hugoniot relations, presented in Reference 3. 

A quasi-one-dimensional characteristic line was then drawn from 
this new shock location through time to an old inner point. Flow properties 
were obtained, and the characteristic line was computed. The information 
obtained was used in the compatibility equation given in Reference 3, and 
flow field variables obtained from the compatibility equation were then 
compared to values obtained from the Rankine-Hugoniot equations. A new 
shock velocity was guessed until these flow field variables converged. 
Thus, an iterative process was established: properties based on a given 

shock velocity were obtained from the Rankine-Hugoniot equations; 
different values were obtained from a characteristics approach and were 
determined independently by a compatibility equation; another shock 
velocity was guessed, and the process was repeated until flow field 
properties at the shock converged. Finally, another time step was taken, 
and the entire process was repeated. 

A similar computaiton was performed for the body. The process 
was simplified, since the body points did not move with time. 

The governing equations utilized for the inner points between the 
shock and body are presented in the form given in Reference 3. 

The spatial derivatives on the right-hand side were evaluated 
using finite differences. MacCormack's predictor-corrector technique was 
utilized in determining the flow field properties after a given time step, 
and details of the procedure can be found in Reference 2. 

Finally, for an upper boundary condition beyond the sonic line, the 
flow is supersonic. That is, any information obtained there will not affect 
properties lying below the sonic line, and upper boundary values may then 



simply be linearly extrapolated 


All flow field properties were therefore computed In this manner 
for each time step taken. 
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1 t l \ A 1 i* j I 
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c 

C Program STUB Jerry Yen 

C Version 5.83.SS 

C 

C This program will utilize the time-dependent technique to constru 

C the shock structure generated by a supersonic blunt body. 

C 

C Limitations: 

C First-order accuracy for spatial derivatives 

C S - D analysis 

C 

C NOTE: M,N = 8,8 

C 

c 

PROGRAM STUB 

COMMON p i , gamma , v i nf , p i nf <. r o i nf , a i nf * k , 

u ( S , 3 ) , v ( S , 8 ) , p < 8 , S ) , ro ( S , 8 ) , sp ( 3 , 2 ) , bp < 8 , 2 ) , b ( S > * 

- dbdy(4) 


Need to input actual numbers for m,n in ALL COMMON blocks 
Also, need m,n for REAL block below 

Rea 1 pi, gamma , v i nf , p inf ,ro inf , a i nf , k , 

. u , v , p , r o , sp , bp , b , dbdy , 

. epss 9 epsb , epst , del , del y , del t , wtO ( 4 ) , theta ( 4 > , 

. cape (3,8) % capp ( 8 , 8 ) ^capr(9 ? S) , capb ( S , 8 ) , psi (8,3) , 

. 2 © ta ( 4 ) , del ( 4 ) , pause ? del zet 

I nt eqer rn , n , step , i cho i ce , j k , 1 1 

prints , 7 Please input number of grid points <n,m): 7 

r tad * , n , m 

print*- , 7 F lease input number of time steps: 7 

read* , step 

call I MOAT ( m , n , ep ss , ep sb , epst ) 
c a 11 ST AR 7 ( in , n , d e 3 >: , d e 1 y , d e 1 1 ) 

Do 1 u j = 1 « m 
1 0 w 1 0 ( j ) = 0 . 

Do 50 jl: = 1, step 

call A T SMi ( w 1 0 , d e 1 1 , ep s s , t h e t a ) 
ca 1 1 AT BODY ( m « n , del y , de 1 1 , epsb ) 

call CAP t m , n , theta «de 1 x , wtO , cape , capp , capr , capb , psi # 

- 2 e t a , d e 1 ) 

call AT I NTR ( cape , capp , capr , capb , ps i , de 1 y , de 1 zet , del , del t . 
call E<NDRY(m,n) 
p a use " s t e d / j k 

if (pause . oe . 0 . 25 . and - pause. 1 1 .0 . 5 ) then 

P r i n t * , 7 Number of t i me step s : 7 , j k 

print*- , 7 Do you want to continue 7 7 
print* , 7 C 1 1 Yes C 81 No 7 
read* , i cho ice 

if ( l cho i ce . eq . 1 ) then ORIGINAL PAGE IS 

call 0P(m.n,jk) OF POOR QUALITY 

else 

goto 60 
end if 

el se i f < pause . ge . 0 . 5 . and . pause . 1 t . 0 . 75 ) then 
print* , Number of time steps: 7 , j k 

print* , Do you want to continue? 7 
print* , 7 til Yes TP3 No 7 

read* , i cho ice 
if (ichoice.eq.l) then 
ca 1 1 OP ( m , n , 1 1 ) 



g .1 sei f ( pause -ge.O. 75 - and .pause. 1 1 . 0 . 9 ) then 
print* , Number of time steps: * « j k 

print* , " Do you want to continue? ? 
print* , ’ Cl 3 Ves CEJ No •' 
read* »ichoice 
if ( ichoice.eq. 1 ) then 
call OF' ( m j. n ? 1 1 ) 
else 

goto 60 ORIGINAL 

end l f op poor 

end i f 

50 continue 

60 Print* , ? END OF PROGRAM ’ 

Print* , ’ Number of time steps = ? , j k 
call OFF IN (m,n, step) 
s t o p 
end 


page is 

Quality 


c 


C— 

C 


C= 

c 


c 

c 

c 

c 

c 


SUBROUTINE INDAT 


Subroutine INDAT ( m , n , e p s s » e p s h « e p s t ) 

COMMON p i % gamma , v i nf ? p i nf , r o i nf , a i nf 

u ( B , S > , v < S « 9 ) , p ( 8 , S > j, r o ( B , 8 ) , sp ( 8 , E ) , bp ( S , 2 ) , b i 8 ) • 
dbdy i h ) 
integer rru n 

Rea 1 pi, gamma , v i nf , p i nf ? ro i nf ? a i nf , k , u , v .. p , r o * sp , bp , b , 
epss ? epsb , epst , ar g , dbdy 
open ( 3 * f l 1 e= • shock . dat ■ ) 


Upstream conditions and number of grid points are inputted 

(m refers to the total number of y points, n refers to the 
total number of >: points); k csfers to a:-: i symmetric geometry 

i ead ( , * ) vi nt , p i nf , r o i nf , g a mm a , k , ep ss , epsb , epst 


C 

C 

C 

C 

C 

c 


20 


C 

C 

C: 


c 

c 

c 


30 
h 0 


a rq — a a m m a * p i n f / r o i n f 
a i n f = S D R T ( a r g ) 

j i oft-ji s t -O the y location! inputted are the ;; locations for 
the n i von j ( y location) ; thus , the arrays contain a set of 2 
numbr.i : t.»e and y location for the body and shock in 

t h e p h y s i c a 1 p 1 a n s 

Do 20 j -- 1 , m 

rend ( 3 . * ) ( bp \ i , i ) , i - 1 , 2 ) , ( sp ( j , i ) , i = 1 , 2 > , dbdy ( j ) 

continue 

Initial guesses for u* p, and ro are inputted; i and j refer 
to points in the physical plane equidistant from each other; 
thus « the arrays contain one number: the value of the property 

a t the loc a tion ;sy (i,j) 

Do m *0 j = 1 , m 
Do 30 i = 1 , n 

r e a d ( 3 • * ) u ( i , i ) , v ( j , i ) , p ( j , i ) , r o ( j , i ) 
cont i nue 
cont i nue 
c 1. o soO } 
r e t u r n 
e n d 



n n n 


Subroutine START ( m , n , de? 1 x , de 1 y . de 1 1 ) 

COMMON p i . gamma . v i nf , p i nf ? r o i nf «. a i nf . k 9 

LI ( 8 , S ) , v ( S . 8 ) , p ( a ? 8 ) , ro ( 8 , 8 ) ? sp ( 8 . E ) . bp ( 8 . 2 ) , b ( S ) , 
dhdy ( 4 ) 

Rea 1 pi? gamma . v i nf . p i nf P roi nf . a i nf . k « u «. v . p , ro ? sp . bp . b . 
d b d y ? 

del x < m ) . delv ( m ) , del t . del t x ? del ty ? arg . a 
Integer n.m 
delt = 1000.0 
Do 20 j = 1 . m 

de 1 x ( j > = ( 1 . O/FLOAT ( n— 1 ) ) *■ ( sp ( j . 1 ) -bp ( j . 1 > ) 
de ] y ( j ) = bp ( m . 2 ) /FLOAT (m-1 ) 
b < j ) b p ( j t 2 ) 

Do JO i = It n 

pd,j) = p(i,j)/pinf 0 RIGINAI da 

ro ( i , j ) = ro ( i , j ) /ro i nf 0F '!? L PA 

u > i > j ) = u ( i « j ) * SQRT ( ro i nf /p i nf ) QL/, 

v* ( i t j ) = v ( i t j ) * SORT ( ro i nf /p i nf ) 
arg - gammat-p ii o> /ro ( i , j ) 
a = SORT (arg) 

d e 1 t x = d e 1 x ( j ) / ( u < j . i ) + a ) 

deity = d e 1 y ( j ) / ( v ( j «. i ) + a ) 
if (deltx.lt. deity. and . d e 1 1 x • 1 1 . de 1 1 ) then 
delt = delt m/1. 5 

e 1 sei f ( del ty • 1 t • de 1 t x • and .deity. It. del t ) then 
delt = del ty/1 .5 

end i t 
c o n t i nu e 


C O iTTTnl 

return 

end 


SUBROUTINE CAP 


Sub r on 1 1 ne CAP \ ,r> , n . theta . de 1 x . w . cape . capo . capr % capb . ps i . 
z f ? t a i d e? 1 ) 

COMl 101 J png a mm a • v i nf , p i nf . r o i nf . a i nf , k . 

u (8«S) . v ( 3 . 3 ) t p ( S . S ) . r o ( 3 . 3 ) . sp (9.2) . bp (3.2) . b ( 8 ) % 
dbdy >. 4 ) 

Rea 1 p i . q aroma . v i nf . p i nf . r o i nf . a i. nf . k . 
u . v . p . ro • sp . tp . b ? dbdy , 

t h g? t a \ m ) .del x m ) . w ( m ) . c ap c ( m . n ) . c ap p < m . n ) .capr (m y n) . 
capb(ni.n) »psi ( m . n ) . z e t a ( m . n ) .del ( m ) 

I n t eg or m . m 1 . n . n 1 

ml -m-1 
n 1 =n— 1 

Do 20 j “ 1 . ro 

Do 10 i ~ 1 . r* 

del ( j ) ~ sp i j * 1 ) - bp ( j . 1 ) 
z e t a ( j , i ) = ( i - 1 ) *d e 1 x < j ) 

cape (j.i) = ( zeta( j . i ) — 1 . 0 ) * d b d v ( j ) - (:e ta(j.i)/ 

TAN ( theta ( j )) ) 
c.app(j.i) = AL.OG ( p ( j . i ) ) 
capr (j.i) = AL.DG (rc<( j . i ) ) 

capb (j.i) = (u ( j , i ) -w ( j ) *2eta ( j , i > +v ( j . i ) *capc (j.i)), 
(sp( j.i) -bp (j.i)) 



n n n n 


r t u r n 
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C 


c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


c 

c 

c 

c 

c 

c 


This is good for converting inner point parameter values into 
values at a certain location 

Subrou t i ne OF* ( 


10 

EG 


>: a nd y 1 o c a t ions 

bp < j * 2 ) = y location (at 
b p ( j * 1 ) . = m 1 o c a t i o n f o r 
n 1 - n— I 

x locati c< n o f p a r tic u 1 ar 
Do 20 j = 1 * m 
y < j ) = 

d x = ( 1 .0/ (FLOAT < nl ) ) 
Do 10 i = 1 * n 
x < j * i ) = ( n- i ) * dx + 
cont inue 


j > 

g i v en y 1 o c a t i o n 
g r i ri p o i n t at y ( 

* ( sp ( j ? 1 ) - bp ( j 

bp < j s* 1 ) 


< j > 

= b p ( j * 2 ) ) : 
, 1 ) ) 


Thus , the value at that x and y location is: 
u ( j j ;: > 

The shock’s x value at a given j/y location is: 
ep < j * 1 > 
r e t; u r n 


c 

C SUBROUTINE ATSHK 

C = = === = =: = = ^ = ^ = = ^ = = = = = =. = = = 


c 


Su b r o u t i ne AT EHK ( w fc 0 * d e 1 1 ? ep s s ? t h © t a ) 

COMMON p i « gamma * v i nf * p i nf * r o i nf « a i nf * k * 

. u ( 8 nB > n v < 8 * 8 ) «i p ( 8 * 3 > * r o ( S • B ) * sp (3»2) * fa p ( 8 * 2 ) «* 

. d b d y ( ) 

Rea 1 pi* gamma * v i nf * p i nf * r o i nf * a i nf * f . u * v * p ? r o * sp * bp * dbdy * 
C 

C * ■* * * * t J 0 T E : The n & : : t two 1 i n © s r e q u i r e #r* ' s a ctual value 
C 


c 

c 

f 

c 


spnew\E*2) * wquess ( -i ) ,ushkO(^) *pshkO(H) .roshkO(H) * 
vsh I 'J \ <+ ) * 

the t a ( m > * a r n ? w 1 0 ( m ) jdeltjepssjus h k * p sh k * r c« sh ! : * v sh k * 
d ©d t * p t ap < 2 ) * p t a ( 2 ) * d © 1 y * d e 1 :: * i ncr y * incrz * 

' ' ypos 1 « ypos2 * zpos 1 * zposE * uy 1 * uy2 * vy 1 * vy2 * py 1 , py2 * roy 1 , 

r oy2 * ual * val * ua * va * pa * roa * aa > p tap 1 ( 2 ) « comp x % comp y * 
rhsal * rhsaE * rhsa3 * rhsa * rhshk 1 * rhshk 2 - r hshk;3 * rhsh k ? rhs ? 

• ucmpat*dif'f 

Integer j?jj*jjj*ii * m * n * m 1 * n 1 * i c h o i c e 

The shod.: point is moved* according to wtO, and at the same y 
location; also found is theta* the angle between the tangent 
to the shock and the x-axis 


Do 1<.> j — ■ 1 9 m 

soncw ( i • 1 > = wtO(f) * delt + so ( i“* 1 ) 



n n n n n n nn 


10 


jt (jij.eq.l) then 
theta ( 1 ) = pi/2, 
g o to 5 

el seif (jjj.eq.m) then 
3 j j = m - 1 
end i f 
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arg - t spnew ( j j j + 1 , 2 ) - spnew ( j-1 , 2 ) ) / < so new ( j j j + 1 9 1 ) 
- spnew ( j-1 , 1 ) ) 
theta ( j ) = ATAN(arg) 

cal 1 EGSHK(wtO( j ) , theta ( j ) «,ushl<0( j ) ? pshkO( j ) , roshkO( j 
vshkO ( j ) ) 

conti nue 


The Ran!, me Huconiot equations are used to determine lu p, 
ro on the shock 

wciuess(j) = 1.0 
Do 20 j = 1 < .n 

ca 1 1 EQSHK ( wguess ( j ) «. theta ( j ) «. ushk ? pshk , roshk , vshk ) 
art} = g aroma -*pshk /roshk 
a s h k =- s Q R T (a r g ) 


A characteristic is drawn from _the new shock point to a po; 
and u . p, and ro are interpolated at fi(tO) 


1 1 0 


120 


120 


dedt = ushk - ashk 

p tap ( 1 ) = spnew(j.l) + ded t* del t*SIN< theta v j ) ) 
p t a p (2 ) = ?pn?w( j,2) - decJ t *de 3 t *C0S ( that a ( j ) ) 
pta<l) ~ (ptap(l) - bp ( j , 1 ))/( spnew ( j , 1 ) - bp ( j , 1 > ) 
p t a ( 2 > = Ao5 ( p tap < 2 ) ) 
m .1 = m - 1 
n 1 = n - i 

i i = 1 
j i ~ 1 

dely - bp (m, 2) /float (ml ) 

* dele -- 1 . 0/ f 1 o a 1; ( n J ) 
i nc r y = p ta ( 2 ) / < de 1 y* j j > 
if ii nc r y . ge . 3 . ) then 
j j - j j + 1 
no to 120 

e .1 'ri e 

i nr: r :: = p t a ( 1 ) / ( d e? 1 z i i ) 

i f < ] nc r s.ge. 1 . ) then 
i i = i i + 1 
cd to 130 

end 3 f 

end i f 

y p o si = d e i v ^ i j 
y p o s2 = d c? 1 y * < j j ~ 1 ) 
zposl = dele * ii 
spos2 = dele * \ii-l) 

J J " J J + 1 

i i ' = ' l i + 1 

if < j j . g t - m . o r . i i . g t . n ) t h en 

print* , 7 THERE IS A PROBLEM IN ATSHi 7 
print* « Please begin with a new w 7 
go to 1150 
end i f 

call I f 4 TER P ^ z p o s 1 • z p o s2 ? p t a ( 1 ) . u ( j j . i i ) , u ( j j * i \ — 1 ) , uy 1 ) 
c a 1 1 I NTERP ( epos 1 » zpo s2 « p t a i 1 ) , v ( j j < i i ) , v i j j « i i - 1 ) , vy 1 ) 
call 1 f 11 ERP ( epos 1 . cpos2 , p t a ( 1 ) , p ( j j , i i ) , □ ( j j « 1 1 - 1 > . py 1 ) 
ca] l I NTERP ( zdosI , cpos2.pt a ( 1 ) ? rc« ( j j , i i ) , ro v j j , i i - 1 ) , r oy 1 ) 
cal 1 I NTERP ( zposl . zposE . p to ( 1 ) • li ( j j - 1 , i i ) ,u( j j-1 , i i — 1 ) 
ca ] 1 : NT FRP < zoos 1 , rnos2 « o t a M ).v( i i - 1 .ii \ . w ( , , - i ; ; _ i \ v ' - 


> 9 


and 


nt A 



n n n n n n d n 


0 . • o s 1 « ypos2 , p t;a ( 2 ) ? uy 1 ? uy2 « u a 1 ) 

call ! : L '* - ' '• p o 5 1 - y p d s2 ? p t a ( 2 ) * vy 1 ? vy 2 ?va 1 ) 

ua = Li - 1 * SINN theta < j ) ) - val * COS < the t a ( j > ) 
va = ual *■ CDS ( theta ( j ) ) + val * SIN ( the t a ( j ) ) 
ca .1 1 I NT ERF’ ( ypos 1 ? ypc*s2 ? p t a ( 2 > ? py 1 ? py2 ? pa > 

call INTERP ( yposl ? ypos2 ? p ta < 2 ) ? roy 1 ? roy2 ? roa ) *Gi) 74L P/\Q?r jg 

a rg = gamma * pa/roa ^ POOR QUALITv 

a a « SGRT(arq) lY 

A new characteristic slope is computed? and A locations are 
evaluated and compared? and iterated until convergence 

dedt =s (dedt + (ua - aa) )/2. 

ptapl(l) = spnew(j?l) + ded t*de 1 t *S I N ( theta ( j ) ) 
ptapl(2) = spnewC j?2) - de?dt*de?l t *COS< theta < j ) ) 
comp = ABS ( p t ap 1 i 1 ) - p tap ( 1 > ) 
co/hpy " ABS(ptapl(2> - ptap(2)) 

1 f \ c omp ; : , g t . ep s s - o r • c omp y . g t . ep ss ) then 

p t ap \ 1 ) = p t ap 1(1) 

p t ap ( 2 > = p t ap 1(2) 

q o to 11 0 

& 1 se- 
pt a ( 1 > = i p tap 1(1) -bp < j ? 1 ) ) / ( sp ( j ? 1 ) -bp ( j ? 1 ) ) 
p t a •* 2 ) - AB3 ( p t ap 1 ( 2 ) ) 

end i f 


The right side of the compatibility equation is computed at 
A a n d a t t n e n ew s h o c I p o i n t 


C 

c- 

c 


1 000 


1 1 00 


rhsa 1 - (va^aa /gamma ) •* ( ( py 1 — py 2 ) /de ] y ) *S I N ( thet a ( j ) ) 
rh s a 2 = aa-k ( u tv 1 *C 0 S ( tine t a ( j ) ) +vy 1 *S IN ( thet a ( j ) ) ) - ( uy 2 * 

COS ( the t a < j ) )+vy2*SIN( theta < j ) ) ) ) *SIN( theta ( j ) ) /dely 
r h sa3 — va * ( v uy 1 * S IN ( theta ( j ) ) — vy 1 *C0S ( theta ( j ) ) ) — ( u y 2 -*• 

S I M \ th et a ♦. j ) > ~vy2* CDS ( theta ( j > ) ) >*SIN< theta < j > > /dely 
r h * s a ~ “ 1 „ v i r ; r t. i : ol+ r h s a 2 — r h s a 3 ) 
j i j ~ i 

3 f ( j j j . p q ,1 ) the n 
rhshk - 0 . 
qo to 50’/ 

els i f ( j i j . a q . m ) •; h (.? n 

j U - m “ 1 

end j 1 

r hi s h l 1 - ( a » ** v s •*. » 0 ( j ) / q a mm a ) * ( psh k 0 ( j j j + 1 ) -psnk 0 v j j j — I ) ) 

S i I* I ( t h o r c . ■ j ) ) / d e 1 y 

r hi s h k c J — a •: h I ♦*.'/£' r* ‘ : *.) ( j j j -h 1 ) — v s hi 1 : 0 ( j j j — 1 ) ) *S I N v t h e t a ( j ) ) /del v 
rhs h ! 2 = vr-h ? - usr.hO i j j j +• 1 ) -u? : h I 0 ( j j j-1 ) ) *SIN ( theta ( j > ) / del , 
r hi s i 1 1 = -1.0^ r hi . . r, i l + hi s hi \ 2- r h := hi 13) 
rht- ~ ( j . t.i ( rriSf. “ rhshik) 


NGUESS is evaluated and a new guess is inputted if desired 


Lt c in p a t “ ( a = hi h ^ t- L 0 G ( p s hi / p i n f ) - a a ■* A LOG (pa/pi n f ) ) / g a m m c 
- - + ua - v rhs ) *■ delt 

d i f f = ABS ( uernoa t - ushk ) 
wr i te ( * ? 1 000 ) itsh I. . uempat ? d i f f 

format!//? U from Rank i ne-Hugoniot = ? ?E10-6?/? 

U from compatibility = ! ’?E10.6?//? 

? hi f er t?nc e = • ? E 1 0 - 6 ? / / ) 

p v i n t* , ? • 

print* ? 7 CONTINUE ITERATION (t ID YES; C2T l-JO ) ? p 

read* . i choice 

if ( ichoicc ,cq. 1 ) then 

wr ) if (*■•!! 01 '* w guess ( j ) 

i or nut ( / » No { -.t recent guess for w = ? *E10.6?/> 


n n n 





call EQSHK ( wguess ( j ) y theta < j > ? ushk ? pshk y roshk ? vshk ) 
u ( j n ) = vshk*COS < theta ( j ) ) + ushk*SIN ( theta ( j ) ) 
v(j,n> = vshk*SIN( theta < j ) ) - ushk*COS < theta ( j > > 
p(jyn) = pshk 
ro( jjn) = roshk 
wtO ( j ) = wguess \ j ) 
sp < j , 1 > = spnew ( j n 1 ) 
sp(jy2) = spnew( j ? 2) 


end i f 
cent i nu fi- 
re turn 
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C SUBROUTINE ATBODY 

C 

Sub r Du 1 i no ATE CD Y v m y n ? del y y del t «. epsb * ph i ,ca pp « c apr ?capv) 
COMMON p i , g anuria y vinf y p i nf « ro i nf < a i nf , k «. 

- u(S?3j • v ( 3 , S ) ? p < 8 y 0 ) y r o ( 8 y 8 ) « sp ( 3 » 2 ) «. bp ( 3 «. 2 ) , 

dbdy ( 4 ) 

Re a 1 p i y g am.Tia * v inf ? p i nf yroinf ? ai nf > \ v.p »ro ysr.p ybp yd b d y * 

. phi ( rn ) ,capp(m) , capr ( m ) . cap v ( in ) y del y * dpdy y drdy * dr dt y 

. d vd v y v 1 « d i f f , de 1 1 , epsb 

I n t e g e r m « n « j j j 
Do 1 0 j = 1 • m 
J J = j 

1 1" i j j j . eq • 1 ) then 
phi ( J ) = pi/2, 
q o i: o 5 

c-? 1 s o ] f < j j j . eq . m ) t h en 
j j j = m - 1 
end 3. f 

ar q = ( sp < i j j + 1 y 2 > - sp(j-l ? 2))/(sp(jjj+l,l> 

- sp ( j-1 y 1 ) ) 
ph 3 ( j > - AThN ( arq) 

5 cap n C j ) = ALOG(p < j 9 1 ) ) 

c a p r v i '> = ALOG (rc( j , 1 ) ) 
cap v i j ) = v ( j y 1 ) /S I M ( ph i ( j ) ) 

10 continue 

J = 1 is evaluated already 

Du 20 j ” 2y IT, 1 

p ( j y 1 > — p N » 1 ) - ( r c* i nf * v i nf#v i nf ) + (roi nf * v i nf * v i nf * 

. COS i ph : v j ) ) **2 ) 

dpdy = ( capp i + 1 ) -capp ( j-1 > ) / ( de 1 y*2 . 0 ) 

dr dy - ( cap r \ i + 1 ) -capr < j-1 > ) / < de 1 y*2 . 0 ) 

10 drdt = ( 1 . O ' gamma )■*( dpdy— gamma *drdy ) *v ( j ? 1 ) 

d vd y “ ( c c.p \ v j + 1 ) -c ap v < j-l ) > / < d e 1 y *2 . 0 ) 
c^pr(j) - capr(j) + drdt*delt 
ro ( j « 1 ) = \ 2 . 71 S23 1 828 ) **capr ( j ) 

d vd t =* — v i i * 1 > *d vd y — ( p ( j y 1 ) / r o ( j ? 1 ) ) *dp d y *S I M ( ph i ( j ) ) 
vl = (capvi j) + dvd t*del t ) *SIN < ph i ( j ) ) 
d i f f = ABS ( v i j y 1 ) —v 1 ) 
if ( d iff . gt . epsb > then 
v ( j y 1 ) = v 1 

qo to 15 

else 


n n n 


V ( j , i - V 1 

u < j , 1 ) = V ( j ; 1 ) *TAN ( ph I ( J ) ) 

end i f 
£0 continue 
return 
end 


ORIGINAL, PAGE IS 

OF POOR QUALITY 



C 


SUBROUTINE ATINTR 


ATINTR c a 1 c li 1 a t © s i n t e r i o r p o i n t properties, 
c ci p r { i « j ) * Lt ( i , j ) , v( i j j ) j p s i ( i , j ) 


Note that input values need to be transferred to subroutine; 
also., variable Cm, n) needs to have m and n transferred; 
don't forget that (m ?n) refers only to INTERIOR points 

Subrout i ne ATINTR < cape , cap p , capr , espb ,psi«,dery,delzet;.del , 
d e 1 1 ) 

COMMON p i , g a, Tina , v i nf , p i nf ,rcinf,ai nf , k , 

u ( 3 , b ) « v ( 8 * B ) , p ( 8 , S ) ? r Tt ( 3 1 S ) ,sp(S,E) ,bp(8,2) , 
d b d y ( h ) 

Rea 1 pi, gamma « v i nf , p i nf ? r o i nf , a i nf , k , u , v , p , r o , sp * bp , db d v , 
c a p c ( rri , n ) « c app ( m ? n ) jeapr (m»n) jcapbimin) , p s i ( m , n ) «• 
ric 1 r © t , del y % de 1 ( »p ) , de 1 1 , dr d z , d r d y , dud z , d udy > d vd z , d v d y 
dps ] d: * dp si dy 9 dpdz , dpriy , 


C* **** * *N0 TE : The ne:: t 3 line require nr" s actual value 

C 


dud t (5,8) ,dvdt (8,8) ? dpsid t (8 , 8 ) ,drdt<8,3) ,rbar (8, 
ubar (3,8), vbar <8,8 ) , ps ibar (8,8) , pbar (8,8) , caprnw ( 
u new i £ , 8 ) , vnew (8,3) , p s i new (8,3) , 

d u d t- 1 , ri v d 1 1 , d s d 1 1 , d r d t; 1 , ri r d t a v , d u d t a v , d v d t a v , d sd t a v 


In 

t e Q sr 

n , rii , 

n 1 , m .1 

nl 

- n 

- 1 


d e ] 

L zet ^ 

= 1 . 0 

F L C A T ( n .1 ) 

m 1 

= m ■ 

~ 1 


Do 

2':) j 

= 1 , 

(li 1 

Do 10 

i “ 

5 , nl 


Fo ) all i 1 1 1 c? r i o r v i , j ) , spatial dsr i va 1 1 v & s a re calc u 1 a t ed 


cal 1 FUD (capr ( j , i ) , capr ( j + 1 , i ) ,del zet , drdz ) 
cal] FUu ( capr < j , i ) , c ap r < j n i + 1 ) , d e 1 > , d rdy ) 
call F WC \ li ( j , i ) , u ( j h* 1 , i ) , de 1 z e t , d u d z ) 
call FUE ; u ( j , i ) , li < j , i + 1 ) , del y , duriy ) 
call FUD ( v ( j , i ) , v ( j + 1 , i ) , de 1 z© t , dvdz > 
ca] 1 FI JD ( v ( j , i ) , v ( j , i + 1 ) , d e 1 y , d vd y ) 

* call FUE ( psi \ j , i ) , ps i ( j+ 1 , i ) , de 1 zet , dp si oz ) 
call FUE .psi ( j , i ) , p s i ( j ? i + 1 ) , d t? 1 v , dp s i d y ) 
call F UD v p ( j,i ) ,p( j + 1 , i ) ,delzet,dpdzeta) 
call FUD v p ( j , i ) , p ( j , i + 1 ) , del y 9 dpdy ) 

c a 1 1 EQMQT ( j , i , c a p c ( j , i ) , c a pp < j , i ) , c apr ( j , i ) , c apb \ j , i ) % 
p s l ( j , i ) , d e 1 ( j ) , dr d z , d ud z , d vd z , dpd z , d ud y , d vd y , d r d y , 
dpdy • dps l dz , dps i dy , dud t ( j , i ) ,dvdt( j,i ) ,dpsidt( i,i ), 
d r d t < j , i ) ) 

rbar( in) = canr(.i,i)’ + (drdt<j,i) * delt) 


ro od 



n n n 


vljcii » J i i 7 — V l _ ■ . i . i . J . I . i . U t_‘ 1 l - > 

psibar ( j , i ) -- psi < j « i ( dpsid t ( j . i ) * del t ) 

pbar ( j . i ) = psibar ( j . i ) + gamma-fc-rbar ( j . i ) 

10 continue 

80 continue 

J = 1 is first analysed; the rest are analysed accordingly 
Do 15 i = Si nl 

call REAR < rbar ( 1 s i ) j rbar < 1 * i — 1 > * del zet j drds ) OR/GfNAL PAGE IS 

call REAR ( rbar ( 1 . i ) * rbar C 2 . i ) . del y . drdy ) OF POOR OliAi itv 

call REAR < ubar < 1 . i > . ubar ( 1 . i — 1 ) * de 1 zet . dud s ) ** * ' 

call REAR ( ubar ( 1 . i ) . ubar ( S . i ) «. de 1 y * dudy ) 

call REAR ( vbar ( 1 • i ) • vbar ( 1 «. i -1 ) ? del set , d vds ) 

call REAR ( vbar ( 1 . i ) . vbar ( 8 « i ) •» de 1 y «» dvdy ) 

call REAR ( psibar ( i . i ) . psi bar ( 1 . i -1 ) . de 1 set . dpsids ) 

c a I 1 REAR ( psibar < 1 . i ) . p s i bar ( 8 . i ) . d e 1 y * dp s i d y ) 

call REAR ( pbar ( 1 . i ) , pbar ( 1 * i - 1 ) . de 1 set - . dpdseta ) 

call REAR (pbar ( 1 . i ) . pb ar ( 8 . i ) . d s 1 y . d p d y ) 

cal 1 EQMOT (Inin cape < 1 . i ) * capp ( 1 * i ) . capr 1 . i ) . capb ( 1 . i > . 

• psi < 1 . i ) .del ( 1 ) . dr dz ? dudz . dvri s . dpd s . dudy . d vdy . drdy . 

. dpd y . dp) si ds . dps i dy . dud 1 1 .dvdtl , dsd 1 1 , drd 1 1 ) 

drdtav = 0-5 * (drdt(i.i) + drdtl) 
dudtav = 0.5 * (dudt(lpi) + dudtl) 
dvdtav = 0.5’* (dvdt(l.i) + dvdtl) 
dsdtav = 0.5 * (dpsidt(l.i) + dsdtl) 
caprnw( 1 . j ) = capr ( 1 . i ) + (drdtav *• del t ) 
unew(l,i) = u(l.i> + (dudtav * delt) 
vnew ( 1 . i ) = v ( 1 . i ) + ( d vd tav *- delt) 
psinewilji) = psi(lj»i> + (dsdtav * delt) 

15 continue 



call REAR ( rbar ( j . i ). rbar ( j . i - 1 ) , del set . dr ds ) 
c a J 3 REAR ( rbar ( j . i ) , r b a r ( j — 1 . i. ) . d e 1 y . ri r d y > 
cal 1 REAR ( ub a r ( j . i ) . ubar ( j . i — 1 ) . cf-e lset .duds ) 
c a 3 1 REAR < ub ar ( j . i ) . ub ar ( j - J . i ) . de 1 y . dudy ) 
call REAR ( vbar ( j . i ) . vbar ( j „ i - 1 ) .delsetjovds) 
ca 1 1 REAR ( vbar ( j . i ) . vb ar ( j - 1 . i ) .del y . d vdy ) 
col I REAR (psibar ( j . i ) . ps i bar ( j . i -1 ) .del set . dps i ds ) 
c a 3 3 R E A R (psi b a r ( j . i ) . psibar ( j — 1 . l ) .do 1 y . dpsidy ) 
c a 1 1 REAR ( pbar ( j . i ) . pbar ( j . i - 1 ) . de 1 zet . dpd seta > 
c a 3. 1 REAR ( pbar ( i . i ) . pbar ( j - 1 . i ) . del y * dpd > ) 
call EO k AjT \ j , i , cape ( j . i ) .capp ( j . i > .capr ( j . i ) . capb ( j . i ' . 
psi ( -,i . i ) .dal v j ) . dr ds . duds . d vds . dpd s . dudy . d vdy . d r dy . 
dpd v .dpsids . dps i dy . dud 1 1 . dvd 1 1 . d sd 1 1 . drd 1 1 ) 


drdtav 

- 0.5 


( d r d t ( j . i ) 

+ 

dr dt 1 ) 

dud t av 

= 0 . 5 

* 

(dudt; ( j , i ) 


d u d 1 1 ) 

d vd t a v' 

- 0 . 5 


( dvd t ( , i ) 


dvdtl ) 

d s d t a v 

= 0.5 

-jpt* 

( d p s l d t ( j . i 

) 

+ dsdtl 


c a p r rwt < j • i > = cap r ( j . i ) (drdtav * delt ) 
Ltnc-M ( j , l > = u(j.i) + (dudtav * delt) 

vr.ew ( j . i ) ~ v ( j . i ) + ( dvdtav * delt) 

- - psinc'W<j.i) = psi(j.i) + (dsdtav * delt) 
30 continue 

40 continue 

Do 105 j - l.nU 
Dr- 90 i = 2 . nl 

u ( j . i ) ~ utifw ( j . i ) 

v ( j . i ) = vr.ew ( j . i ) 

r d ( j . i ) = < 2 . 718831 828 ) ** \ capr n;; ( j . i ) ) 
arg = psincv;(j.i) + gamma^capr nw ( j . i ) 
p(j.i) = (2,71 988 1 228 ) **arg 
90 Cent i nue 

1 05 C:«nt i nue 


enci 


ORIGINAL PAGE IS 
OF POOR QUALITY 

c 

C SUBROUTINE BNDRY 

C 

Subrout i ne BNDRY (m,n) 

COMMON p i q aijuria <• v i nf 9 p i n f « r o i n f ? a i rrf n k ? 

. u ( S i S > * v ( Q •> 8 ) • p ( 8 ? S ) 9 r o < S n 8 ) * sp ( S « 2 ) « bp < 8 • 2 ) ? 

. Tdbdy<4) 

Real pi * gamma ? v i nf* * p .i nf « r o i rrf * a i nf ? k « u ? v « p «. r o ? sp % b p «. dbd y 
J n t c? u S i (Ti 9 n 





Do 10 i = 

1 * n 




u ( m « i ) 

= ( ( u ( m- 1 !« i ) —u ( m-2 * l ) )/ui ro-2 * i ) ) * 



a 


u ( m- 1 « i ) + u v m- 1 - i ) 




v v in O ) 

= ( ( v ( m- Ini) -V ( m-2 «. i ) ) / v ( m-2 9 i ) ) * 

_ 


a 


v v m— 1 9 i ) + v ( m— 1 «. i ) 




p v m * i ) 

= ( ( p ( fTi 1 ? i ) -p ( m-2 % i ) ) /p < m-2 n i ) ) * 



a 


ptm-l si ) + p(m-l«i) 




r o *' m - i 1 

1 = ( v ro ( m~l , i ) -to < m-2 <. i ) ) /ro ( m-2 n i > ) * 



a 


ro ( rn- Ini) + ro ( m- 1 * i ) 



1 0 

Co nt i nue 



r e t u r n 


end 


C 

C SUBROUTINE FWO 

C 

C FUD generates the forward spatial derivative at ( i ? j ) : 

C dvard • : i j 

C 

Subroutine FWD(vari j *varil j « di-: ? dvard;: i j ) 

I ’i rv a 1 v* r\r 1 ) <* van 1 j n d j * d v a r d y * i j 

d v a r d r. i j “ \ v a r i 1 j — v a r i j ) / d m 

r e t u r n 
end 


C 

cr SUBROUTINE REAR 

C = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = ======= 

C - . 

C REAR generates the rearward spatial derivative at i?js 

C dvd::ijl 

C 

Subr out i ne REAR *, var i j , var i ml j « d.x « d vd x i j 1 > 

Rea 1 var i j < var i rr? 1 j « d>i <i d vd x i j 1 

d v d ■: i j 1 = ( v a r i j — v a r i m 1 j > / d x 
rr-tur n 

GTid 

c 


n n n n n n 9 


sL'r-r.L'u i I InE riTTERF 7 


INTERP interpolates for a value answer 

Subroutine INTERP < r owl , rowE , mi d , co 1 1 , co IS » answer > 

Real rowl . rowS.mid < co 1 1 j co IS » answer 

answer = ( <mid-row2> * (col l-col2 > > / <rowl-row2> + colE 

return 
end' 


OF 


Po °» 


SUBROUTINE EOMOT 

C, 

C EQI10T computes values for time derivatives of the 

q equ a t i o n s o f hi o t i o n at ( i ? j ) J 

C drdt, dud t . dvdt, dpsidt 

C 

Subroutine EOMOT ( j - i , cape , capp .capr , capb , psi . de 1 .drdi.dud:, 

dvd:: , updz , dudy ,dvdy , dr dy ,dpdy , dpsidz , dpsidy , dud t , dvd t , 
. dsdt jdrcJt ) 

COMMON p i • aniiiiiid n v i nf i p i nf ?r o i nf nf 

u < 3 ,3 ) , V ( 8 , 8 ) , p < 8 , a ) , ro < S , B ) , sp ( S . 2 ) , bp ( 8 . 2 > , 

. dfcdyt4> 

Real p i , gamma , v inf , p inf , r o i nf , ai nf , k , u , v , p , r o , sp , bp . dbdy , 

. cape .capp , capr , capb , psi > del , dr dz , drdy , dudz ? dudy , dvdz , 

dvdy , dpsidz , dp sidy , dpdz , dpdy , dud t , dvdt , dps id t , drdt , y 
I n taper 1 1 j « On n 
y - 1.0 

drdt = - ( ( c apb * drdz ) + (dudz/del) + 

. < (cape/de 1 ) * dvdz) + dvdy + <v< j , i)*drdy> + 

. < k* ( v ( j • i ) /y > ) ) 

dud t = -< ( capb*dudz ) + (v( jui >*dudy >+< (p< j > i >/<ro< j « i )^rlel ) ) 

. * d p d z e t a ) ) 

dvd t =- ( (capb^dvd z ) + ( v < j < i ) *dvdy ) + < i < p < j * i > *capc ) / ( ro ( j > i ) * 
del)) * dpdzeta) + ( ( p ( j , i ) /ro ( j > i ) > * dpdy)) 
dpsidt = -((capb *• dpsidz) + < v ( j ? i ) * dpsidy)) 
r g? t u r n 
end 


C 
C= 

c 
c= 
c 
c 
c 

Subroutine EQSHK(w, theta «ushk f pshU ?roshk ? vshk ) 

COMMON p i . gamma « v i nf ? p i nf * r o i nf * a i nf ? k « 

, „ u(SsS) ,v 1*8,8) !.p<S, 3 ) ,ro(8jB) ,sp< 8 , 2 ) •bp(S ! .S) * 

. -dbdy ( 4 ) 

Rea 1 pi ..gamma • vinf * p inf ?roi nf ? a inf ? k ? u ? v ? p ? ro ? sp 5 bp ?dbdy ? 
w theta «unel «. ushk > pshk ? roshk •» vshk 

unci = ( gamma- 1 . )*< < vinf-w ) *SI N ( the ta ) >**2+E.*ainf*ainf 
us hk - unel /\ igamina+1 - >* < vinf-w > *SIN ( theta) )+w*SIN< theta) 
pc-.h I- = (£.* ( <vinf-w)*SIN< theta) ) ( gamma- 1 . ) )/<gamma + l . ) 

roshk = ( <gamma+i . > *pshk+ < gamma- 1 . ) )/( < gamma* 1 . ) + 

. ( gamma-1 . ) *pshk ) 

vzhk “ vinf * COS (theta) 
ret ur n 
end 


SUBROUTINE EOSHK 
Rani i ne— Hugo mot rel at i ons 


c 

c 

c 

c 


SUBROUTINE OP 


OP outputs values in the physical plane 


OFd!, NAL PAs E is 
0F P °0 R quautv 


1 000 
1 020 

1040 


/ / / • 


, / / ) 
,/ ) 


Sub r o u t i ne OP ( m * n 9 1 1 ) 

COMMON p i , g amnia , v i nf , p i nf > r o i nf * a i nf , k , 

u < 3 , 8 ) , v ( 3 , 8 ) , p ( Q 8 ) , ro ( 8 , B ) , sp ( 3 * 2 ) « bp v 3 < 2 > » 
dbdy < 4 ) 

Rea 1 pin gamma v i nf «. p i nf « r o i nf « a i nf * k , u , v , p « r o , sp ? bp ? dbdy • 

C 

C******** NOTE! the next line requires actual m^n values 
C 

. m ( 8 n S ) j y ( 8 ) n d >: 

Integer t t 

□pen i 4 ,file=- shock . out ? ) 
wr i t e ( 4 % 1 000 ) 1 1 

format ( / / , ’ T 1 ME STEP NUMBER : / > 

wr i t e v 4 < 1 02 0 > 

format (///* 31 OUTPUT FROM STUB 

’ PHYSICAL LOCATIONS OF GRID POINTS: 
wr i te ( 4 « 1 040 ) 

format (jX , •" IX J V 

DO 20 j - in ITl 

v ( i ) = bp C j «. 2 ) 

dr: = (1 .0/ (FLOAT (n-1 > > >*<sp< j , 1 )-bp < j * 1 > ' 

Do 10 i » 1 «. n 

;■: ( j , i ) = (n-i) * dr: + bp<j>l) 
wr i t & < 4 n 1 050 ) i * x ? j « y 
format ( 5 X’ , 14, 2 X nF r 8.4n2X ? 1 4 , 2 X , F 8 • 4 ) 
cont i nus 
cont i mil? 
wr i t e ( 4 • 1 030 ) 

format (//,' PHYSICAL LOCATIONS OF SHOCK AND BODY POINTS: % 
//,-• dJotir that j: coordinates are presented for a given y > ’ • 
5X , ' Y X locations: SHOCK BODY ’ »/> 

Do 4 0 j *= 1 « »Ti 

wr j te ( 4 « 1 080 ) y ( j ) • 5 p ( j ? 1 ) ? bp ( j ? 1 ) 
form nt ( 5X • F8 . 4 * J. 6X * F3 - 4 * 3 X , F8 . 4 ) 
c o n 1 1 nu c 

wr i te \ 4 , 11 00 ) 
f or mat ( / / « 

.PROPERTIES: ? , / ■ > 

Do 30 j = 1 « iT> 

wr i 1 1 ? \ 4 , 1 1 5 0 ) y ( i ) 

fc r ma t i / / , ’ Y LDCAT I ON : LFS.4, / / ) 

ur i tc i 4^ 11 75 i 

format •: 5*. , 7 X U V p R0 ' ' n 

r.'c- 50 i - 1 , n 

'wr i t D i 4 » 1 200 ) x ( j » i ) ? u (jji)jv(j«i);p(j«i)jroCjji) 
f o r ivu*. t ' 5 X * FS • 4 ? 2 X i 4 < F 8 - 4 « £ X ) ) 
cont :i vine 
cont i ni.ie 
C 1 DEB V 4 ) 
return 
end 


1 050 
10 

20 

1030 


1 080 
40 

1 1 00 


1 1 50 
1175 

* 

1200 


PHYSICAL LOCATIONS AN! VALUES OF F^ulJ 


C SUBROUTINE CPF IN 

C 

c 
c 


OP output'.: va 1 lies in the physical plane 


c 

c******** 

c 


1 000 
1 020 


1040 


1050 

10 

£0 

1 060 


i oao 

40 


LI ( 8 , B ) , V < S , 8 > , p < a , 3 ) , r o ( S , 8 ) , sp < 8 , £ ) , bp ( 8 , 2 ) . 
dbdy(4) 

Rea 1 pi .gamma , v i nf . p i nf . ro i rrf . a l nf , k , u , v . p . ro , sp . bp . dbdy . 


NOTE: the ne:;t line requires actual m?n values 


:•: ( 8 , 8 ) ? y ( 8 ) ,d>: 

-Integer tt 

JDpen ( £ . f i 1 e= ’ shock .fin 7 ) 
wr i t e ( E . 1 000 ) 1 1 

fo mat ( / / j ' T I ME STEP NUMBER : 
wr i te ( £ , 1 OHO ) 


0F poor 


pag* 

Qvalit 


fs 


, 14.// ) 


format <///,' OUTPUT FROM STUB ’,///, 

7 PHYSICAL LOCATIONS OF GRID POINTS: './/) 
wr i te ( 2 -1 1 040 ) 


.. -format <5X, ’ I X J Y ’,/> 

Do 20 j = 1 • m 

y \ j ) = bp ( j *2 ) 

ch: =*■ < 1 .0/ (FLOAT <n-l )> >*<sp < j , 1 >-bp< j, 1 > > 

Do 10 l = 1 n n 

i ) = ( n—i ) * d >: + bp ( j ? 1 ) 
wr i t f? ( 2 9 1 050 > i ? « ? j ? y 
.format <5X, 14 ,HX ,F8.4 ,2X * I4,2X ,FB.4) 
cent i nue 


conti nue 
wr i te ( 2 9 1 060 ) 

format <//,’ PHYSICAL LOCATIONS OF SHOCK AND BODY POINTS: 

/ / %’ (Note that :: coordinates are presented for a given y) 9 % 
5X , 11 Y X locations: SHOCK BODY ; ? / ) 

Do 40 j = 1 , m 

wr i te < 2 9 1 080 ) y(j);isp(j ? l) 9 bp(j 9 l) 
format ( 5X , FS . 4 , 1 6X , F8 . 4 « 3X , F3 . 4 ) 
co'nt i nue 


wr i t q ( 2 9 1 1 00 ) 

1100 format ( // < PHYSICAL LOCATIONS AND VALUES OF FLOW 

.PROPERTIES: ' ,// ) 

Do 60 j =■• 1 « m 

wr i te ( 2 1 1 1 50 ) y < j ) 

1150 format (// , ? Y LOCATION: ’ ,FB.4,//> 

write*: 2 , 1 175) 

1 1 75 forma t ( 5X • ’ X U V P RO > , / > 

Do 50 i = 1 * n 

wr i t e v 2 « 1 £• jO ) >: < j ? i ) ? u ( j 9 i ) 9 v ( j 9 i ) 9 p ( j 9 i ) ? r o ( j 9 i ) 

1 200 f o r ma t < 5 X < FS . 4 , 2 X 9 4 ( FS . 4 9 2X ) ) 

50 continue 


APPENDIX A3: SEAGULL 


Description of the Program SEAGULL 


The program SEAGULL is a code designed for the analysis of 
two-dimensional or axisymetric supersonic inviscid flow of an ideal gas. 
All the discontinuities are treated explicitly. The limitation of the 
program is that the axial component of the mach number should remain 
supersonic during the execution. Although originally designed for internal 
flows, SEAGULL can be used for external flow due to optional input of the 
lower or upper wall geometry. The geometry for the upper wall and lower 
wall are assigned in subroutines GEOMB and GEOMC. Descriptions of the 
input and outputs are given on the following page. The code could also be 
used to merge several ducts into one analysis. 

For the purposes of the propulsion group, the program is used to 
calculate the flow properties along the forebody of the hypersonic drone. 
The floMBwpwrafwe^s obtained *• from this program will be used for the 
boundary layer code EDDVBL to obtain the shape of the boundary layer. 
Once the boundary layer geometry is known it could" be added to the 
forebody geometry for a better estimation of the flow properties along the 
forebody. 


The input parameters 
Steps between output 

Number of sharp corners on upper wall and lower wall 
(10 maximum each)* their locations and angles in degrees 
Type of the flow (2D or axisymetric) 

Mesh points 
Geomatry discription 
Constant or variable Gamma 
Mach number 

Distance to end the run 

Free stream Pressure* Temperature* Gamma* Inclination in 
degrees, Mach number, Gas constant 
Lower and upper wall temperatures 


PylByt description 

Forces and Moments on the wall and on the R-Z plane 

Viscous and Inviscid forces 

Heat transfer from upper and lower wall 

The pressure* temperature* velocity componant in the R and Z 
direction, the mach number, the streamline inclination in 
degrees, total temperature, for each mash point and its R 
location. 
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Verification: Seagull 


For the five degree cone, the downstream values output from SEAGULL 
at Mach 10 upstream conditions are: 

M - 8.271 

P/P. - 3.060 
RHO/RHO. -2.139 


From [22], for a five degree cone with no angle of incidence to the 
freestream and M^-10, we have: 

M « - 9.925 

P/P. - 2.95 

RHO/RHO.- 1.8033 

A comparison of these values illustrates a limitation of the code; 
specifically, SEAGULL'S downstream values are somewhat in error. Also, 
SEAGULL shows that the flowfield downstream of the shock decreases in 
velocity with increasing distance from the nose. These effects are caused 

by the duct flow modelling that SEAGULL uses. This effect is minimal, 
however, and the properties on the wall change by less than 5% over the 
entire forebody if the walls of the duct are kept far apart to minimize the 
compression caused by the internal modeling of SEAGULL. This introduces 
some error into the downstream flow properties. This error is unavoidable 
due to the necessity to model multi-ramp geometries, which was not 
possible by merely using conical shock tables because the ramps are not 
subject to uniform freestream conditions. 
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PRO GRAM EDDYBL 

FDDYBl in a f u llv compressible, two-dimensional or axisymmetrir 

computer proq-'am for the calculation of boundary layer properties. The 

program accounts for mass flux at the body surface, heat flux. 

pressure and temperature gradients, and both transverse and 

longitudinal bodv curvature. The program has the capability to start 

* - 

from a Stagnation ooint and proceed through transition j r , t: 

turbulence. Th* equations used by the program cannot account, 

however, fo- the- effects of shock-boundary layer interaction or 
sepa rated flow. If the program encounters separation, execution 
terminates at that point, although the data is saved for all 
calculations prior tc that point. Also, the program does not account 
ior na ' : lficatfrion or if's effects in any wav. 

9 

7h ® ‘ a • u 1 a 1 1 o •*. of the boundary lave. properties cur. be . r .* c .V : 

• ••. w-j.s. utiuw s* :.t,r «.r. ' s t ws-equa t ion .-odcl or it’s mult isoala 

( \»'Vt 

mode. . T i ef . i i ■ ] . In short. the two equation model uses conservation 
of hkv, and momentum i conjunction with an equation for the turbulent 

• i nq energy, and another equation for the turbulent dissipation 

ij 

r « » -• . C ■ e f . *t . d. 51. These last two equations account for the creation 

'If-?. t i* ur t ) On of 1 u* h\j 1 ,c (7 . 

The multi scale model differs from the two equation mode! ;r, 
•" v|llil the mult i -scale model also uses the ennse. i cn 

C:f mass and momentum, as well as the same turbulent m,*inrj ene.g-. 
equation', the multiscale model uses a slightly different equation for 
the turbulent dissipation rate to “predict the effects of plane strain 
and shear on homogeneous turbulence." fr P f^. p.q]. The larqe *.. 

d.f fcarencu between thy two models lies i n the wav that they treat the 


R •’* j r i o 1 <i ' ^ i Stress t pnsor , 


which accounts for vortices and turbulence. 
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The two-equal; ion mode: Iuiuds all of the vortices, or '’eddies", ir.to 
on., term. The mu literal*? model Jumps the eddies into "large" and 
small " eddip... I , p.*T]. This scaling of the turbulence proves 

^ ‘ ^ ‘ m > u bl u.irrjl..; fh&n t wo - equ i I i ort mod* 1 1 . Clue app ] ications use t h-~ 
j J l i ..c a 1 « f-m/.u-J 1 'j^ivr-iy, and sr further d i sc us?, i o n will be 

limited to the multi scale model , References 1 through c offer a much 
more cun.f.lete dc-sc- lotion ot the two models used. 

As t-reqo'.is)y mentioned. fTDDVFL is a fullv compressible boundary 
larer code. In addition to the other capabilities previously 

mentioned. EDPVBL has the canability tc start form user - supp 1 i ed 

hour, bar. layer oo-?r t i es . That is, if the user knows what the- 

properties ot the boundary lave* are at a certain point in the 
f ! owt i«ld. then c omr.utat io *. c*n begin at that point and there is no 
"~ e ' : to .t i< 9 Io r a : f u 1 a t i ons for points further upstream. Fo<*lhet more 
this allows restarting FDDYHi. from previous runs. In c-ur application, 
wc- ‘"11 "SC EOr>'«BL t.-. calculate boundary lsye.- properties to a point 
who' e the i r c ■ stream conditions change. Ue then restart computation 
from the old boundary lave-r properties with the new -free stream 
n .id i t l c»ns> . Jo tiiis wav wr- proceed > r. a pi ere wise manner along the 
h" dv i r. r< i.i-i 1 i i> . 

The.-u am .-racy input p v am • e. s f or FDnvSl . The free sti ..-am 
yf'io-ltv, static temr-c-r a l.-jr ■ • .nd pr e-isun? ml-, be tilhe- t-; input 

dl.-ectl.vv or mav fcp . a 1 o. 1 a t M hyrhe program itself. The progr -nr. will 
Cdlrulali.- the free, stream properties hae.r-g „ pon the far upstream total 
pressure, total t^pivat.jrrj. mn:h number „nd shock wave angle. •' EDDVEl 
“-• e ” U,r ' P9nl i -v-Hun.-jni o 1. fV. I a f ions to calculate the free stream 

if.r rl« ngi i«‘* t n of I h? ). 
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The properties of the free stream gas itself are also defined by 
t.he user. The ratio of specific heats* the specific gas constant. and 
both the 1 am i n,tr sod turbulent Prandtl numbers must be defined by the 
The turbulent Prandtl number has been determined experimental lv 
^ t»e 0.9 for. most flows, but communicat ior- with the author indicator 
that more accurate results are obtained if the turbulent Prandtl 
nun.be* is set to 8/9. The laminar Prandtl number is a function of 
viscosity, which is in turn a function of temperature. It's value be 
calculated by the user and input foi each free stream condi t. i on. 
Another free stream proper ty that is user-supplied is the free stream 
turbulence intensity. This value varies from 0.05V. for well designed 
wind tunnels to 1 00V for turbulent jets. Cref.l^. p.P57D. 

The properties of the surface itself are described by three- 

jWMHIUHIIIIII lil HI ill ■ ' 

separate data files. The first data file contains N ar c 1 eng th-pr essurr- 
pairs that define the static pressure at the body’s surface. (The user 
determines the number of arc 1 eng th-presssure pairs. N). Specifying the 
P r ^ssui’ e in this manner al lows the user to impose pressure gradients 
upon the boundary layer. EDDYHL uses a cubic spline fit to determine 
pressures at intermediate points. 

The second file defines the thermal properties of the surface. 
Again. f*> set?, of arc length, temperature and pressure values are 
specified. Note that only temperature or heat flux need be specified. 

If temperature is specified, heat flux will be calculated. If heat. 
■Hu* is .specified, then temperature will be calculated. 

The last data file specifies the body geometry of the surface and 
the rate of mass addition or removal at the wall. N sets of five 
numbers define the mass flux and body geometry. Fpur of the numbers 
describe the arc length, axial distance, the radius from the 


center 
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line, and the local curvature for a given arc length. The fifth number 
specifies the mass flux at that station. 

F i na 1 1 v . within the program itself the user defines the angle 
thf#* the "free stream mates with the center line of the body. Also, the 
roughnG^ss height Reynolds's number is defined to allow for different 
surface roughnesses. 

The grid resolution for the finite difference method is also 
determined by the user. Initial and final arc lenghts, maximum number 
of stream wise steps and the initial step size are all defined bv 
the user. The step sir. e between stations will be changed bv the 
program itself as needed to maintain accuracy and conserve time. 
F ><rC u t i c n w i 1 ] be terminated when either the specified arc length has 
been the? program has reached the* maximum number of stream, 

wise steps. Also, the user can specify the number of grid points 
normal to the surface. More points will increase the resolution but 
will also increase computing time. The number of points normal to the 
surface remains constant. Also, the user must choose the model to be 
used. The closure coefficients for the models mav also, be changed, but 
1 1 e rtutho f has shown the default values to perform well when compared 
to experimental results, [ref. n ,\s i . Also, the velocity profile of the 
boundary layer obevs the Sutherland Viscosity Law, althouah a power 
law m-'iv l.*o invoked bv the user. 

Finally. the user must define boundary layer* properties if 
starting from his own pr ede t ^r m i ned c o nd i t i ons . The coefficient of 
' - » f r_i r t l u* ■ • the* boundary 1 t h i c I ness , shape factor, and 

' ■; • ■« l -JiS-. d i; » • »» t • i* . 1 L : must a 1 1 he ) i . J , 

the u'.iPr in addition to the parameter*:, mentioned 
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The program outputs boundary layer properties for each stream 
wise station, t see attached sheet and ref X, op. 53-2^ ]. The .output 

listing is very complete and easilv lends itself to comparison with 
f-Kpei imental ^ 5 ta. Add i t i ona 1 1 y , the program output-3 a boundary 1 aver 
profile at the last station and at any intermediate station requested 
bv the 1 user. It alsG creates an output file that can be used to 
restart from the last station's properties. 

The program has several basic limitations. First, it is limited 
to laminar or turbulent, non-rarif ied boundary layer analysis- A: 
previously mentioned, the program cannot account for separation o** 
shock * boundary layer interaction. Also, free stream values remain 
constant throughout a given run, and if the free stream is indeed 

variable then the program must be run in a series o' stream wis.= 

» 

sections, each section maintaining constant free stream c ond i t j c n=> . 
Finally, the program is limited to t wc-ri imens i ona 1 analysis, although 
the effects of stream wise curvature are accounted for in t u r 
a x i s v mmr? tr ic case. 



Vprifiration: EDDYBL 


Program EDDYBL is a product of DCW Industries and has been tested by 
the author for compressible flows ranging as high as Mach 5 and including 
the effects of cooled walls and adverse pressure gradients. [19, p.11-13]. 
Accuracy of the computed values is shown to be within 8% of experimental 
results. 

Despite verification of the program's accuracy, there are a number of 
factors affecting the program's application to hypersonic boundary layers. 
Among these are surface roughness and freestream turbulence. Another 
possibly important aspect which was not considered in this application 
were real gas effects, which become important at high Mach numbers. The 
groups original intention was to compare the results of this program with 
those of another boundary layer program which accounted for real gas 
effects. This comparison was never accomplished. However, results in- 
literature for high Mach number flows and communication with researchers 
at NASA indicate that the freestream conditions under consideration are 
not sufficiently severe to warrant consideration of real gas effects. [16]. 
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APPENDIX: BDDYBL OUTPUT PARAMETERS OF KOcH QUAU7Y 


The following quantities are printed in the integral parameter 
portion of Program BDDYBL output. 


Name 

Svmbol/Equation 

Definition 

Dimensions 

BETA 

F=(2?/U e )dU e /dC 

Pressure gradient parameter 

None 

CFE 

c fe = 2^»/ Pe^e^ 

Skin friction based on edge density 

None 

CFW 

c fw =2 ' c w/ Pw^e^ 

Skin friction based on wall density 

None 

D/T 

H=S*/9 

Shape factor 

None 

DELTA 

S 

Boundary- layer thickness 

ft 

DLTAST 

S* 

Displacement thickness 

ft 

DPEDS 

d(p e /p.U. 2 )/dJ 

Dimensionless pressure gradient 

None 

DTBDX 

dT e /dT 

Transformed edge temperature gradient 

None 

DUEDX 

dU e /dT 

Transformed edge velocity gradient 

None 

ECOUNT 


Number of points where k < 0 

None 

EMAX 

JtH )ny/ 

Maximum turbulent energy parumeter 

None 

ES 

s 

Dimensionless surface dissipation rate 

None 

ETAEDG 

N e 

Mesh point number at B.L. edge 

None 

HD 

fi=q w /(T w -T aw ) 

Heat transfer coefficient 1 

Btu/ft 2 /sec/°R 

IEDGE 

N 

Total number of mesh points in B.L. 

None 

ITRO 


Number of iterations for convergence 

None 

M 

m 

Streamwise step number 

None 

ME 

Me 

Edge Mach number 

None 

MUB 

Me 

Edge molecular viscosity 

slug/ft/sec 

NSKIP 

- 

Number of points below u = USTOP 

None 

NSTE 

fi/p e U e C p 

Stanton number based on edge density 

None 

NSTW 

fi/Pw^eCp 

Stanton number based on wall density 

None 

NUE 

Pr^hs/MgCp 

Nusselt number based on edge viscosity 

None 

NUW 

Pr L hs/M W C p 

Nusselt number based on wall viscosity 

None 

PE 

Pe 

Edge pressure 

lb/ft 2 

OSD 

<Jw 

Surface heat flux 

Btu/ft 2 /sec 

RE 

Pe 

Bdge density 

slug/ft 2 

REQBLT 

R*S> = PeU e $*/ Me 

Displacement thickness Reynolds number 

None 

RES 

Re,=Pe u e s/ Me 

Arc length Reynolds number 

None 

RETHET 

Re^p^e/ng 

Momentum thickness Reynolds number 

None 

RFC TOR 

r 

Recovery factor 

None 

RMI 

r o 

Body radius 

ft 

RVWALD 

Pw v w 

Surface mass flux 

slug/ft 2 /sec 

S 

8 

Arc length 

ft 

TAUD 

Tw 

Surface shear stress 

lb/ft 2 

TE 

T e 

Edge temperature 

°R 

THETA 

e 

Momentum thickness 

ft 

UE 


Edge velocity 

ft/sec 

UTAU 

u r 

Friction velocity 

ft/sec 


Name 

Symbol/Bouation 

Definition 

Dimensions 

WCOUNT 

* = /Pe u ePe r o 2 ^ da 

Number of points where w < 0 

None 

XI 

Transformed streamwise coordinate 

None 

Z 

' z 

Axial distance 

ft 


The following quantities are printed in the profiles portion of 
Program BDDYBL output. 


Name 

Symbol /Equation 

Definition 

Dimensions 

i 

i 

Mesh point number 

None 

y/delta 

y/s 

Dimensionless normal distance 

None 

u/Ue 

u/U e 

Dimensionless velocity 

None 

yplus 

y + =u T y/\i w 

Sublayer distance coordinate 

None 

uplus 

u*=u/u r 

Sublayer velocity coordinate 

None 

k/Ue**2 

k/U e 2 

Dimensionless turbulent kinetic energy 

None 

omega 

v e c^n e 2 

Dimensionless dissipation rate 

None 

eps/mu 

PT/m 

Dimensionless eddy viscosity 

None 

L/delta 

<7TK/(uS) 

Dimensionless turbulent length scale 

None 

uv/tauw 

<-pu'v’>/x^ 

Dimensionless Reynolds shear stress 

None 

T/Te 


■ ■ Temperature ratio 

None 
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A computer program entitled "Design of Supersonic Inlets by 
a Computer Program Incorporating the Method of Characteristics" 
(Ref. 11J, dubbed "Inlet" by the propulsion group, was obtained 
to aid in the design of the inlet of the scram jet. Inlet is a 
sophisticated program written in 1968 by Bernhard H. Anderson of 
the NASA Lewis Research Center in Cleveland, Ohio. Inlet is 
based on steady, inviscid, irrotational , two-dimensional or 
axisymmetric flow. Its main solution scheme utilizes the method 
of characteristics, which was originally developed in 1929 but 
only practically implemented after the advent of modern comput- 
ers. The method of characteristics as used in Inlet has its root 
in two-dimensional potential flow theory. Not only is the poten- 
tial function, phi, included (phi is actually incorporated into 
the method of characteristics equations), but the stream func- 
tion, psi , is also included for determination of the diffuser's 
physical boundaries by numerical mass flux integration. The 
oblique shock relations are also extensively utilized to deter- 
mine the flowfield; hence continuity, momentum, and energy 
equations are also considered in the calculations. Additional 
equations incorporated which yield design consideration data are 
mass flow spillage and additive drag equations. The numerical 
techniques used by Inlet to solve these equations simultaneously 
include finite difference techniques, Runge-Kutta numerical inte- 
gration, and iteration. 

Inlet is particularly useful in designing scramjet inlets 
since it a-llows the internal boundaries of the cowl or center body 
walls along with freestream conditions to be prescribed. This 


results in the calculation of the internal cowl and centerbody 
wall geometries required for the input flow conditions. A useful 
application of Inlet is to calculate the internal cowl and 
centerbody boundaries from given upstream flow conditions and 
desired downstream conditions at the entrance of the diffuser. 
This is done by specifying forebody ramp geometry and freestream 
flow conditions along with desired shock pattern at the cowl lip 
and internal to the diffuser. Inlet then calculates the geometry 
required to cancel or reflect the shock at either the cowl or 
centerbody boundary by using the characteristic equations. Con- 
currently, Inlet determines the placement of the boundary by 
considering mass flux integration and performing iterations until 
the mass flux obtained by the placement of the boundary equals 
the initial mass flux entering the diffuser. In this manner, the 
internal boundaries are constructed and printed out as part of 
the output data. The boundaries can then be used to study the 
off-design conditions for the scramjet inlet just constructed. 

Some important phenomenon to be studied in the next two 
quarters include the boundary layer, shock wave-boundary layer 
interaction, and chemical dissociation in the boundary layer. 
Because inlet does not account for these phenomenon, a comparison 
of the results obtained from Inlet at high freestream Mach 
numbers (M>5) with experimental data in that regime will be 
essential in justifying the use of Inlet for the design of the 
forebody and diffuser geometry. Such a comparison performed in 
the Inlet description report (Ref II) revealed good Inlet results 
for Mach numbers less then 2.5. However, comparisons at M = 3.85 


ft 


begin to show deviation between the theoretical and experimental 
data, probably due to the increased boundary layer effect at this 
higher Mach number. Since boundary layer effects and shock wave- 
boundary layer interaction become increasingly more significant 
at Mach numbers greater than 5, the determination of the extent 
of this deviation for Mach numbers at the des i'g n region is essen- 
tial if Inlet is to be useful for the final design of the fore- 
body and diffuser. Nonetheless, should deviations prove to be 
significant, Inlet should at least yield some qualitative guide- 
lines to consider in the final design of .the scramjet inlet. 


9T 


INPUTS 


Freestream Conditions 

AMO Free-stream Mach number 

GA M Ratio of specific heats (Use 1.4 for ideal gas) 


Ramp Geometry 

THETC Initial ramp angle, radians 

BETAE Estimation of ramp shock angle 

ALPHA Deflection angle of second ramp (measured from tangent 

of first ramp) 


Cowl Specifications 


COWLA 

THETAL 


SPILL 


Initial interior cowl angle, radians . 

Location of cowl lip relative to spike tip measured 
from centerline of ramp in radians (Note: THETAL 

is used only when the cowl contour is unspecified. 
Under this condition the cowl lip radius is 1. 


Set SPILL = 
shock wave. 
THETAL. 


0 to locate the cowl at along initial 
Set SPILL = 1.0 to locate the cowl at 


XCOWL 


YCOWL 


x-coordinate of cowl lip position (Used when you want 
to place the cowl lip somewhere besides at the location 
where the ramp shock wave would intersect it. Use when 
inlet geometry is known and off — design performance is 
to be studied. Set XCOWL = 0 to place the cowl lip 
along the ramp shock wave.) 

y-coordinate of cowl lip position (Used when you want 
to place the cowl lip somewhere besides at the location 
where the ramp shock wave would intersect it. Use when 
inlet geometry is known and off— design performance is 
to be studied. Set YCOWL = 0 to place the cowl lip 
along the ramp shock wave.) 


OUTPUTS 


Conditions in the Vicinity of Ramp Oblique Shock 
Program outputs: 

1. position of point for calculations 

2 . Mach number 

3. angle of -‘flow 

4. angle of shock wave 

5. ratio of static pressure behind shock to static pressure in 
front of shock 

6. ratio of total pressure behind shock to total pressure in 
front of shock 

Conditions on Ramp Contour 

1. x-coordinate of point 

» 

2. y-coordinate of point 

3. magnitude of velocity vector 

4 . flow angle 

5. x-direction component of the velocity 

6. y-direction component of the velocity 

7. Mach number 

8. ratio of static pressure to freestream total pressure 

9. ratio of static pressure to total pressure behind the shock 
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SPECIFYING COWL GEOMETRY 

Cowl contour is specified by third-order polynomial of form 

Y = R( 1 , 1 ) + R( 1 , 2) X + R(1 ,3) X 2 + R(1 , 4 ) X 3 

The coefficients are specified in sets from 1 to NR where 
NR is an input variable specifying the number of sets of cowl 
contour coefficients to be read in. Also the starting position 
where the above equation is in effect must also be specified by 
the variable COWL(I) which specifies the X-coordinate to impose 
the specified set of coefficients. 


Centerbodv Specifications 

Similiar to cowl specification, centerbody is express as a third- 
order polynomial of the form 

Y = S(I,1) + S( 1 , 2 ) X + S( 1 , 3) X 2 + S(1 .4) X 3 

The coefficients are specified in sets from 1 to NS where 
NS is an input variable specifying the number of sets of 
centerbody contour coefficients to be read in. Also the starting 
position where the above equation is in effect must also be 

coordinate to impose the specified set of coefficients. 


Calculation of Cowl or Centerbodv Coefficients 


When the Cowl or Centerbody coefficients are to be calculated set 
the corresponding coefficients equal to 0 and instead specify the 
flow conditions within the inlet by specifying the following 
input variables: 


NTHR 

AMT(J) 

THR ( J ) 
ANG( J ) 
NIS(J) 


Number of internal isentropic compression sections. 

Final Mach number specified for the end of the J th 
section, J = 1,NTHR . 

Length of section. 

Final surface angle (radians) of section. 

Number ot surface net points within THR(J). 


qo 


Flowfield Specifications 


START 


DELP 


PRINT 


M 


Initial net point spacing parameter for point loca- 
tions. Flow properties are calculated at each point. 

(A good initial value for start is 0.050) 

Ratio of angular location of last inital data point to 
inital shock angle. The program constructs the initial 
datum line from the point on the surface, BODY (2), to 
a point in the field whose angular location is DELP 
times the initial shock angle. 

Set PRINT = 0 to bypass the output printout for the 
forebody flow field calculations. 

Set PRINT = 1.0 to print flow field calculations. 

Number of points in flowfield calculations. 


Internal Flow 

NTHR Number of internal isentropic conpressions sections. 

ISHK Used for off-design calculations only. 

For on-design calculations set ISHK = 1. 

For* off-design ISHK is index which specifies at which 
centerbody location BODY (ISHK) the shock wave is 
canceled . 

NSHK Number of internal shock waves to be computed. 

Set NSHK = 1 to cancel cowl lip shock. 

NSHK - 1 specifies number of internal shock reflec- 
tions. 


Input 

DELB 

DELU1 

ERROR 

OFF 

NDIM 

NR 

ND 


Parameters 

Initial increment used in calculation of ramp shock 
angle (set DELB = 1.4 x 10 - ) 

Integration increment for conic flow field calculation 
(Set DELU1 = 1.0 x 10"°) 

Convergence parameter (Set ERROR = 1.0 x 10 -4 ) 

Set OFF = 0 for on design calculations. 

Set OFF * 1.0 for off-design calculations. 

Set NDIM - 2 for 2-D case. 

Number of sets of cowl contour coefficients. 

Number of increments in addiditive drag calculations. 


Conditions in the flow field 


1. x-coordinate of point 

2. y-coordinate of point 

3. magnitude of velocity vector 

4. flow angle 

5. x-direction component of the velocity 

6. y-direction component of the velocity 

7. Mach number 

8. ratio of static pressure to freestream total pressure 

9. ratio of static pressure to total pressure behind the shock 

Conditions in the Vicinity of Cowl Oblique Shock 
1« x-coordinate of point 
2 . y - c o 5TTTT — " ■ ' ' - ' 

3. Mach number 

4 . flow angle 

5. shock wave angle 

6. ratio of pressure behind shock to pressure ahead of shock 

7. rato of total pressure behind shock to total pressure 
ahead of shock 

Conditions along Cowl Contour Surface 
1 . x-coordinate 

2. y-coordinate 

3. magnitude of velocity vector 

4. flow angle 

5. x-direction component of velocity 

6. y-direction component of velocity 

7. Mach number 


pressure 


8. ratio of static pressure to freestream total 

9. rati-o of static pressure to total pressure behind shock 

Table of Cowl Coefficients 

The program list cowl coefficients along with the x-coordinate 
position where the coefficients become effective. 

Table of Centerbody coefficients 

The program list centerbody coefficients along with the x- 
coordinate position where the coefficients become effective. 
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APPENDIX CP I PE 


Summary 

CPIPE is a fortran program written in I960 by G.Y. Anderson 
of NASA Langley Research Center. CPIPE is a 1 -dimensional , real 
gas analysis of the combustion of hydrogen in air in a supersonic 
channel (see Figure M). CPIPE was originally intended to perform 
a one step analysis of the entire underbody flowfield, but can he 
modified to only analyze the combustor region. Given the initial 
upstream conditions and other inputs (see Tab 1 e fA ) , CPIPE outputs 
the downstream 1-dimensional flow properties. CPIPE allows for 
several fuel injection points and allows the user to specify the 
fuel Mach number » stagnation pressure and temperature* and angle 
at which it is injected. There arc several other options 
available reqardinq the boundary layer transition point and 
mixing models used (see Options, Tabic M). In addition, CPIPE 
contains some empirical relations (introduced hy Chuck McClinton, 
NASA Langley Research Center) which approximate the likelihood of 
autoiqnition and flameholding. 
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Conservation Equations 


CPJPE uses the mass, momentum, and energy conservation 
equations to perform a one step analysis of the flow properties 
in a supersonic channel. The equations and the control volumes 
they are applied to are shown explicitly below: 
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Te mp era ture 


CPIPE performs a 1— dimensional analysis of the wall heating 
* n a ' supersonic channel. Given the temperature distribution 
along the wall of the channel > CPIPE will calculate the heat 
flux. The equations used to calculate the heat flux are shown 
explicitly below: 
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Abstract 


The aim of this paper is to design and compare active cooling systems that met the 
requirements of a hypersonic drone. The design procedure involved the determination of a 
flight condition to be studied; justification of the use of an active regenerative cooling system; 
and, finally, the optimization of the design of such a system. Several computer codes were 
written to facilitate parametric studies of the design. Two designs are proposed: one utilizing 
the hydrogen fuel as a coolant, and the other utilizing ethylene glycol as the secondary fluid. 
Due to weight considerations, the hydrogen design was determined to be preferable. 
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1 List of Symbols 


a Speed of Sound 
c Specific Heat (Btu/lb R) 
g Acceleration due to Gravity 
h Altitude (Ft.) or Enthalpy(/< 2 /,sec 2 ) 
k Conductivity 

q Total Heat Transfer Rate ( Btu/ft 2 ) 
t Time (Sec) 

A Area (ft 2 ) 

D Drag = 0 .5 Peo V 2 C D A 
L Lift = 0.5 p^V 2 C L A 
M Mach Number = V/a 
P Pressure 
R Radius 
T Temperature 
V Flight Velocity (ft/s) 


SUBSCRIPTS 


B body 
d diameter 
eq equilibrium 
s stagnation 
sh shock value 
SL sea level value 
w wall value 

oo free stream 

6 inviscid edge of the boundary layer 
SUPERSCRIPTS 
mean value 

derivative with respect to time 
SPECIAL SYMBOLS AND CONSTANTS 

g gravitational constant at sea level 

Le Lewis Number = Pr/Sc 

L/D Ratio of lift to drag 

Pr Prandtl number 

Re Reynolds number 

p 3 l sea level density 

a Stefan-Boltzmann constant 
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GREEK SYMBOLS 


e Emissivity 
p Viscosity 


p Density 



2 Introduction 


The phenomenon of aerodynamic heating was first experienced on re-entry vehicles, where 
heating rates on the order of tens of BTU / ft 2 s were encountered for relatively short periods 
of time. In spite of these high heating rates, the short time interval that the vehicles were 
in the critical regime allowed the problem to be managed with passive systems, such as 
ablative surfaces or heat sink structures. 

Air breathing vehicles, generally experiencing less intense heating rates, tend to be in 
the critical regimes for much longer periods of time. The aircraft skin does not possess 
sufficient thermal capacity to act as a heat sink. In addition, requirements for reusability 
and turn-around time rule out ablative surfaces. 

With these considerations in mind, the aim of the Thermal Management Group was to 
design a hypersonic drone cooling system suitable for testing candidate National Aerospace 
Plane (NASP) cooling systems. The hypersonic drone vehicle configuration was the same (See 
as that studied by the 1987-88 UCLA Advanced Aeronautics Design Class. The analysis 
was focused upon the nose region of the drone since the computation of high stagnation 
point heating rates is simpler and better understood than the heating rates i eenring within 
the engine. 

Due to the relatively low heating rates experienced at the drone nose, the use of a heat 
pipe was unwarranted. Instead, the designed system was based on an actively cooled panel 
concept. Although the actual thermal design of the NASP nose region would probably make 
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use of the heat pipe concept, there will exist panels elsewhere that will require an active 
cooling system. The design presented here, then, is simply a test of that concept. It would 
be judicious to point out that this thermal design is by no means meant to be directly 
applicable to the nose of the N ASP. Rather, it is the concept of actively cooled panels that 
is currently under study; the system design just happens to be located at the nose region 
of the hypersonic drone. 

♦ 

The design procedure involved the determination of a flight condition to be studied; 
justification of the use of an active regenerative cooling system; and finally, the analysis 
and optimization of the design of such a system. The aircraft configuration is presented 
first, followed by a detailed description of the design procedure and results. 

3 Determination of Flight Profile and Nose Geometry 

To obtain useful data from flight tests, a suitable flight profile was established. In general, 
the flight condition is governed by two basic parameters: altitude and Mach number. De- 
creasing altitide will increase air density, resulting in increased heating rates and higher 
equilibrium temperatures, particularly at high Mach numbers. Also, with an increase in 
Mach number, specific impulse (ISP) and engine efficiency tend to decrease. Thus, given a 
range of Mach numbers for scramjet operation (5-12) and a rough estimate of the fuel tank 
size (8 feet diameter, 40 feet long), an effort was made to determine a flight condition where 
active cooling was required (due to material limitations), while maximizing the flight time 
for data collection and experiments. 
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A code called Hypersonic Equilibrium Temperature and Heat Flux (HETAQ) [See Ap- 
pendix 2.1] was created to develop a data base of flight conditions. HETAQ varied Mach 
numbers from 5 to 12, altitudes from 70,000 feet to 115,000 feet, and nose radii from 1 to 
6 inches. Lees Approximate Method was used to determine stagnation point heating rates, 
and an energy balance of the convective heat transfer with radiative cooling produced the 
skin equilibrium temperature. A simple routine incorporating the thrust equation was used 
to determine the total flight time available given the geometry of the proposed fuel tank, 
Mach number, and altitude. The code calculated time to thermal equilibrium for aluminum 
and titanium noses at each given flight condition. This was necessary to prove that the 
time required for the craft to reach thermal equilibrium would be less than the cruise flight 
time. 

Once the data base was developed from HETAQ, the group selected a flight condition 
suitable for the entire design team. Specifically, though a higher Mach number was desirable 
for the propulsion group, the resulting shorter flight time proved to be incompatible with the 
Thermal Management Group’s aim of achieving steady-state conditions. From the HETAQ 
data, any speed faster than Mach 9 or 10 would result in the use of too much fuel and 
thus create undesirable heating rates. Therefore, the speed was set at Mach 10, and the 
altitude set at 100,000 feet. In addition, the nose radius was set at one inch, since a smaller 
nose radius results in a higher equilibrium temperature (the temperature that would arise 
without cooling). The higher temperature was desired to justify the use of an active cooling 
system. Indeed, the equilibrium temperature for such a configuration is 4000*R, high enough 
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to require an active cooling system. 


In summary, with a flight profile of Mach 10 at 100,000 feet, and with a nose radius of 
one inch, results from HETAQ determined that the stagnation point and points aft required 
cooling. 

4 Methods for Computing Aerodynamic Heating 

In the hypersonic regime, the aircraft endures the most severe heating rates at the stagnation 
region of the body. A stagnation point is defined as the point on the body where a streamline 
terminates. The stagnation region differs from the rest of the flow field in that the boundary 
layer is at its thinnest. A large change in velocity near the stagnation point means a large 
change in kinetic energy. This lost kinetic energy is converted to thermal energy in the 
form of stagnation point heating. A thin boundary layer means the flow is at free stream 
conditions very close to the body. Heat transfer is thus at a maximum where the boundary 
layer thickness is at a minimum. The larger velocity gradient and the thin boundary layer 
combine to make the stagnation region the area of critical aerodymnamic heating. 

It is important to note that the boundary layer, defined as the region dominated by 
viscous effects, has a finite thickness at the stagnation point. Thus, the flow properties at 
the edge of the boundary layer, required by all of the methods discussed below, are the 
properties behind the bow shock, NOT the stagnation properties. 


Heat transfer rates from flat plates subjected to convection are given by equations of 



the form 


q 3 = /i(Too - T w ), (1) 

where h is the convective heat transfer coefficient. This quantity depends on the 
Reynold’s Number 

The length from the leading edge, L, is zero at the stagnation point. Therefore, instead 
of relying on Re as a parameter, the methods discussed use the velocity gradient at the 
stagnation point. 

A method used extensively is the Fay and Riddell Method. Developed in 1959 at Avco, 
it is accepted as the most accurate formula available, but is also one of the most difficult to 
use (See Ref. 11). The heat transfer rate is given by 


?j = 


1020Pr° V7C p 9 p 9 


(^r)° •W.n.n*. - + (i." - 1)^ 


] (2) 


where 


Le = Lewis Number = 1.4 
m = 0.52 for Equilibrium Boundary Layer Flow 

= 0.63 for Frozen Boundary Layer Flow 

It is apparent that eqn.(2) is quite complex, particularly since terms such as the en- 
thalpies and the velocity gradient are themselves derived from long equations. 
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Lees Approximate Method, presented in 1958 (See Ref. 11), resulted in a reasonably 
accurate (+ or - 12 percent) value for stagnation point heat transfer rate. 

«■ - < 3 > 

For pencil and paper preliminary design, Lees Approximate Method provided the nec- 
essary accuracy with the minimum amount of computation and the least chance for error. 
Lees Approximate Method was also the chosen as the standard method by Hankey (Ref. 

ID- 

Fortunately, the group had access to MINIVER, a software tool that calculates stagna- 
tion point heating rates as well as heating rates at other regions. MINIVER utilizes the 
Fay-Riddell Method among others, depending on the region and flow regime. Between Fay- 
Riddell, Lees and MINIVER, the group was able to verify results quickly. This practice not 
only served to increase the group’s confidence in the results but also to increase the group’s 
understanding of the physical phenomenon. 

So far, stagnation POINT heat transfer rates ( ) have been discussed. Before one 
can speak of heat loads (§~0, the heat transfer rates away from the stagnation point must 
be determined and the frontal area must be known. The heat transfer rate away from the 
stagnation point is given as a function of theta by 

q = q a {COS6fl 2 (4) 
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for 6 up to 60°. This correlation is accurate to within seven percent (See Ref.ll). Since 
five degree half angle nose cones are considered appropriate for hypersonic aircraft (See 
Ref.4), theta varies from 0 to 85 degrees. This range is beyond the sixty degree limit of 
accuracy, but is adequate for a first order approximation. Using the incremental area as 
shown Figure 6, the above equation is integrated over theta to give heat transfer (§~) 
through the nose. Recall that heat transfer rate PER AREA decreases with an increase in 
nose radius. However, total heat transfer increases with an increase in nose radius. This is 
due to the increase in frontal area with larger radii. An optimum would have to be found 
between the size of the nose radius and and the size of the conical portion of the nose. 
Correlations exist that give heat transfer rates for flat plates at angles of attack. Using the 
proper incremental area, these correlations can be integrated to find the total heat transfer 
to the conical portion of the nose. 

5 Mission Selection 

The first step involved the derivation of stagnation point heat transfer rates (? fe) usin § 
Lees Approximate Method for various flight conditions. Heat transfer rates were computed 
in terms of nose radius and wall enthalpy, h w , given by 

h w = 0.23 4(T W ) + (10)~ S (T W ) 2 (5) 

This allowed designers to vary the important parameters of nose radius and wall tem- 
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perature at a given flight condition and determine the heat transfer rate. The following is 
a list of some of the flight conditions examined. Stagnation point heat transfer is given by 

. / 1 \ . 

q> = ~ v r) (6) 

where X and Y are constants FOR THAT FLIGHT CONDITION. 


MACH 

ALTITUDE 

X 

Y 

6 

50000 

74.07 

777.40 


70000 

49.02 

782.70 

8 

50000 

175.75 

1307.02 


70000 

116.19 

1315.11 

9 

100000 

93.24 

1704.68 

10 

50000 

343.30 

1987.30 


70000 

226.95 

2000.45 


A graphical summary of this data is given in Figure 2 for various nose radii and flight 
conditions. In addition, sample calculations can be found in Appendix 1. 

At this point, the worst case must be examined. Provided there exists no active cooling 
scheme and the structure of the vehicle does not act as a heat sink, the final equilibrium 
temperature of the nose must be determined. Even if the heat from aerothermal heating is 
not bled off by a coolant or absorbed by the structure, energy is still lost to the surroundings 
via radiation. There exists some temperature such that the amount of heat coming in due 
to convection is exactly equal to the amount of heat leaving due to radiation (See Figure 
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3). This temperature is called the equilibrium temperature and can be solved for by using 
the following equation. 


q, = to(T* eq - l£) (7) 

ca n be taken to be zero and as a first estimate, e, the emissivity, can be taken to 
be 0.85. For the same flight conditions listed above and a nose radius of four feet, the 
corresponding equilibrium temperatures are given below. 


MACII 

ALTITUDE [ft] 

EQUILIBRIUM TEMP. [R] 

6 

50000 

2225 


70000 

2090 

8 

50000 

3025 


70000 

2810 

10 

50000 

3777 


70000 

3525 


Note that even with such a large and impractical nose radius, the equilibrium tempera- 
tures become too high for a conventional material such as aluminum. Therefore, an active 
cooling system must be designed that can maintain the nose at a reasonable temperature. 
Figure 4 shows equilibrium temperatures at 100,000 feet as a function of Mach number 
and nose radius. Note that at our designated flight condition of Mach 10, 100,000 feet the 
equilibrium wall temperature is about 4500JZ. 
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The next parameter that must be examined is the time required for the nose to reach 
equilibrium, provided no heat is bled off by a coolant or absorbed by the structure. An 
approximate energy balance that gives an estimate of time required to reach equilibrium 
skin temperature is given by 

pCpV {JT W ' C q — Tinitia/) = q s A&t — <J€A(T w ^q Too) At (8) 

Assuming an “eggshell” type structure of uniform thickness, times on the order of sec- 
onds are required to reach equilibrium (See Figure 5). Aluminum, 1/16” thick takes the 
longest to reach steady state. In reality, there would exist some support structure behind 
the shell (spars, ribs, etc.) that would act as a heat sink and lengthen time to equilibrium. 
However, the craft would have to be accelerated to its target flight condition. As it passes 
through the lower Mach numbers, it would experience aerothermal heating. Therefore, 
when the clock “starts”, T{ n i ixa i is already greater than the temperature at launch. Despite 
these complications, it is obvious that a drone at the designated flight condition would easily 
reach thermal equilibrium within the estimated 5.32 minute burn time. 

In summary, 

1) There exists several methods, including the code HETAQ, that allow rapid and ac- 
curate calculation of stagnation point heat transfer rates. 

2) The stagnation point is the area of critical aerodynamic heating. 

3) Correlations for heat transfer rates can be integrated over frontal areas to determine 
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total heat transfer (§“f) for a reasonable first estimate. 

4) Heat transfer rates and heat loads at hypersonic speeds are large enough to warrant 
design of an active cooling system for any vehicle that spends more than even a few seconds 
in these regimes. 

5) Due to propulsion and drag considerations, Mach 10 and 100,000 feet is chosen as 
the target flight condition. 

6) With an estimated flight time of 5.32 minutes at Mach 10, 100,000 feet, the nose of 
the aircraft should reach thermal equilibrium long before the hydrogen fuel runs out. 



6 Analysis and Design of an Active Cooling System 


High heating rates at the nose of the drone dictate that the coolant be routed from the 
tank to the stagnation region and then back along the cone before being routed to the 
scramjet (See Figure 7). The design is broken up into two parts; impingement cooling and 
panel cooling. Impingement cooling takes place in the manifold immediately behind the 
stagnation region where a single jet of coolant impinges on the inside surface of the nose 
and is then routed into numerous tubes. Panel cooling takes place after the coolant enters 
the tubes and is modelled by a fin analysis. 

Both impingement cooling and panel cooling are modelled in the Actively Cooled Struc- 
ture Evaluation System (ACSES), a FORTRAN computer program written by the authors. 

6.1 Impingement Cooling 

A method commonly used to cool the leading edge of a gas-turbine airfoil involves the 
impingement of a jet of coolant against the leading edge’s inside surface. Impingement 
cooling can also be applied to the nose of the hypersonic research vehicle currently under 
study by the UCLA Advanced Aeronautics Design Class. 

In Ref. 13, Livingood presents Nusselt number correlations determined experimentally 
for jet impingement upon a concave, semi-hemispherical surface. Since the nose of the 
hypersonic test vehicle was modelled as a semi-hemisphere, a subroutine program called 
upon by ACSES utilized these correlations in determining the final fluid bulk temperature 


15 



exiting the semi-hemisphere manifold. 


The geometrical parameters of the manifold (see Figure 8), including nozzle diameter, 
d, nose radius, R n , and distance between the nozzle and the nose inner surface Z n , were 
input into the program. These parameters could be varied to optimize the manifold design. 
The fluid bulk temperature entering the manifold from the nozzle (T b>noz ), and the coolant 
mass flow rate through the nozzle, MFR noz , were also input and optimized. Finally, fluid 
properties for two different coolants were used. The program outputted the fluid bulk 
temperature exiting the impingement cooling apparatus, T btOUt . This T b<out value was then 
returned to the program as the inital coolant temperature entering the active cooling tubes 
extending aft of the nose. 

The procedure used to find T btOU t was an iterative one. First, an initial guess was made 
for T bi0ut . Fluid properties were then calculated based on the film temperature Tju m \ 


Tjiim = + 0.5 (T b>noz ■+■ Tout)], (9) 

where T WtTnax was the maximum allowable wall temperature. For the design considered, the 
nose wall was assumed to be at T w<max for the maximum efficiency of the design. T w<mai 
itself was considered a conservative value since it was taken as half of the material’s melting 
temperature. 

The average Nusselt number correlation, N u avei for jet impingement on a concave semi- 
hemisphere, as presented by Livingood (Ref. 13), was used to determine an average heat 
transfer coefficient h ave . The correlation and the resulting average heat transfer coefficient 
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are presented below. 


N u avc 


2.98 i?e°- 585 (|) 110 (^-)- 0007 ; 






( 10 ) 

( 11 ) 


where, 

Re — Reynolds number based on nozzle diameter, d 

D = 2 R n = nose diameter 

k — thermal conductivity of the coolant 
Once the average heat transfer coefficient was obtained, the total heat transfer into the 

coolant was evaluated. The coolant heat transfer was equivalent to the aerodynamic stag- 
nation point heating of the nose minus the radiative cooling. Thermal resistivity through 
the thin-skinned, highly conductive nose was considered negligible in this analysis. Thus, 
the coolant heat transfer was: 


Q coolant ” Qnoae Qradiation 


( 12 ) 


where q nose was a function of the stagnation heat flux as determined by Lees’ approximate 
method of stagnation heating in Hankey’s report (Ref. 11), and outlined in Appendix 1. 
Given the surface area of the nose, the temperature exiting the manifold could be easily 
determined: 


Qcoolant — have A sur f {Tb'Qut - Tb,noz) 


(13) 
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or, 


Tfc.out = , ] -■—+ T t ,n 0 z (14) 

«at;e^*ur/ 

This r(,, ou< was then compared with the initial guess, and the entire procedure was repeated 
until the values for 7i i0ut converged. 


6.2 Panel Cooling 


The coolant tubes were placed equidistant apart on the cone. The distance between the 
tubes L, was therefore only a function of x, the distance from the manifold. This symmetry 
of the geometry and the assumed symmetry of the aerodynamic heat load dictated that only 
one of the coolant panes need be considered. These panels were then divided into slices of 
width and analyzed accordingly. 



The analysis was divided into three modes of heat transfer: convection from the coolant 
to the inside wall of the tube; conduction from the inside wall to the outside wall of the tube; 
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and conduction from the base of the fin through the nose cone wall (which was modeled as 
a fin). Given an initial coolant fluid temperature, the goal of this analysis was to calculate 
the maximum temperature of the panel. 

0.2.1 Fin Analysis Due to symmetry, coolant temperature at a given x location 
was independent of the particular tube. This fact implied a temperature profile shown 
below, with = 0 at the midpoint section of the two adjacent tubes. 


fig b. 



i 

t 


The analysis was then simplified, in that only half of the panel need be modeled by the 
fin. The fin endpoints were the base of the tube and the midline of the two adjacent tubes. 
A detailed derivation follows. 
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X akin = Thickness of skin. 

T 0 = Base of ”fin” temperature, 

XijX? = Station locations: Ax = — X\ 

Performing an energy balance in the y direction: 


■ Assume Top surface adiabatic: Temperature and heat flux constant across width Ax 


qA c \ y - qA c \ y + dy + q»Pdy = 0 


(15) 


A c = cross sectional area = t - Ax 
P = exposed perimeter = Ax 

•Dividing by dy and using the definition of a derivative and Fourier’s Law: q = —k 
the above equation becomes: 

d 2 T 

kt ~dy* + 9 *( x ) = 0 ( 16 ) 

•However, it was assumed the temperature did not vary across the width Ax so an ordinary 
differential equation was obtained (Eqn 17): 


d 2 T 

k, w + ,>(I) = 0 


( 17 ) 
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•Boundary conditions: 


•Integrating twice: 


dT 


= 0 . 


dy 
T = T 0 . 


• y= 2 


• y = o 




2kt 


V + C x y + C 2 


(18) 

(19) 


(20) 


•Applying the boundary conditions, and solving for the constants, Eqn 20 becomes: 


T = 


2kt 


(y 2 - yL ) + T 0 


( 21 ) 


Equation 21 applies anywhere along the fin; however, the temperature of interest was the 


midline temperature of the panel, which was the maximum temperature the panel of skin 


would reach: T l = T crt , 
y= j 


Tcrit = q,(x){—) + To (22) 

T 0 , the temperature at the base of the fin which was the same as the temperature on the 
outside wall of the tube, was calculated from the conduction and convection analysis. 

6.2.2 Convection The total heat that must be absorbed by the coolant fluid was 
equal to twice the amount of heat absorbed over the surface area of the fin, as heat was 
entering the tube from both sides 


Q = L(x)q, Ax 


(23) 
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The flow in the tube was turbulent with a high temperature gradient at the wall. Thus 
a convection correlation given by Incopera and DeWitt(Ref. 9) was used to best fit the 
conditions. 


Nu d = 0.027 Re 4 / 5 Pr 1 / 3 ( ~)° 14 

mu = dynamic viscosity at the mean temperature of the fluid 

mu s = dynamic viscosity at the wall temperature 

R*d = 

valid if: 


0.7 < Pr < 16,700, 
Re d > 10,000, 

4d > 10, 


•The heat transfer coefficient can be calculated from the Nusselt number: 

Nu d k 
h = ~ 


(25) 


•From Newton’s Law of Cooling, the temperature of the wall can be calculated: 


T wa ll = T m + 


Q 


(26) 


The properties were evaluated at T m , the bulk mean temperature of the fluid. T m was 
calculated by using a heat balance. Knowing an initial T m> i n , T(m, out) can be calculated: 

Q 

Tm,out — -Tg + T m>tn (27) 
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Therefore, the temperature of the inside wall was now known and could be used in the 
conduction analysis to calculate the temperature on the outside of the wall of the tube, or 
T 0 , as it was referred to in the fin analysis. 

6.2.3 Conduction 



The equation for conduction heat transfer through a plane was used to calculate the 
outside wall temperature, T a . The heat was assumed to flow through the bottom of the tube 
only. A fin analysis of the circumference of the tube demonstrated a rapid 

decrease in wall temperature after approximately 30 percent of the diameter. Because the 
tube was welded to the nose wall, additional heat transfer would occur. Thus, 50 percent 
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of the surface area was chosen as the effective area of conduction. 


Acond — 


dir Ax 


( 28 ) 


Treating the tube as a flat plate of thickness and thermal conductivity k tu be> the 

heat rate is: 

Q = —A C0nd (T 0 - T wall ) (29) 

%tube 

This equation can be solved for T 0 


T 0 = 


Q^tube 

ktubeAcond 


+ T wa ll 


(30) 


The conduction heat transfer and the convection heat transfer equations were combined 
into a thermal resistance network. Conservation of energy implied the heat rates, Q, used 
in equations (26) and (29) must be equal. 


network 


&- 
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(31) 
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1 


T 0 = Q( 


hAr 


+ 


't'tubc 


k tubc-Acond 


) + Trr 


(33) 


Finally, T a can be calculated from eqn. (33) and Tcrit from eqn. (22). The objective 
of the design was to keep T ct h less than the maximum allowable temperature for the nose 
cone material, T max . 



7 Results/Discussion 


Two design configurations were evaluated. One utilized liquid hydrogen as the coolant 
and the other a 30 percent aqueous glycol solution. The hydrogen configuration was 
composed of supercritical hydrogen flowing directly from a highly pressurized fuel tank to 
the active cooling system design. The hydrogen was then diverted back to the scramjet 
engine. The ethylene glycol design consisted of glycol circulating through a closed loop 
system including a pump, a heat exchanger, and the active cooling system design. 

Both the designs were formulated with three common constraints in mind: the 
critical temperature was to remain less than the maximum allowable material 
temperature; the coolant velocity was to be less than the fluid speed of sound; and the 
pressure loss was to be kept as small as possible. A description of the specific constraints 
applied to each design is presented below. 

The critical wall temperature for both designs was set at 1000 K. This corresponded 
to half the melting temperature of the wall material, molybdenum. This constraint was the 
most critical for both designs. 

In the hydrogen design, the constraint placed upon fluid velocity was the most 
restrictive. The highest value of fluid velocity occurred in the nozzle entering the 
manifold. The fluid temperature in the nozzle was assumed to be equal to 20 K. The speed 
of sound of supercritical hydrogen was then found to be 1675 m/s (Appendix 4). The 
pressure of the hydrogen stored in the fuel tank was assumed to be approximately 3000 psi. 
It was decided to set the maximum allowable pressure to twenty-five percent or 750 psi. 

Similarly, the glycol design constraints were determined. The speed of sound was 
1644 m/s (Ref. 6), and the maximum pressure loss was assumed to be about 700 psi, to limit 
the size of the pump required. An additional constraint placed on the glycol design 
involved the temperature range of the coolant. The freezing point of a 30 percent aqueous 



glycol solution is 257 K; this temperature was chosen as the minimum coolant temperature. 
The maximum coolant temperature was chosen to be 420 K (Ref. 6). 

Once given the design constraints, various parametric studies were performed and 
several observations were made. In order to decrease the pressure loss, large tubes or 
smaller coolant mass flow rates were needed. However, with a smaller mass flow rate, the 
nose temperature increased and more tube branching was required. The larger the tube 
size, the less efficient , 'the heat transfer in the tube. Hence the critical wall temperature 
increased with increasing tube size. 

7.1 Final Designs 

Liquid Hydrogen Co olant Final Design : 

After parametrically varying mass flow rate and tube geometry, the liquid 
hydrogen design converged upon the following configuration (see Appendix 3.1): 
mass flow rate : 3.5 kg/s 
Tube diameter : 4.88 mm 
wall thickness : 2mm 
initial number of tubes : 24 

Branching locations : 0.225 m, 0.875 m, 2.1m, 3.55 m 

Coolant velocities are kept well within the sonic range in both the impingement 
nozzle (806 m/s) and the tubes (281 m/s). System pressure loss is also kept to a minimum of 
about 150 psi. Outlet coolant temperature is 703 K. Therefore, all of the previously 
mentioned constraints were satisfied. 

Increasing tube size would lower the pressure drop even more for the same mass 
flow rate. However, the fewer number of tubes allowed would result in the critical 
temperature being exceeded. Hence, 4.88 mm diameter tubes are the largest standard 
tubing that can be used with out fear of the craft exceeding the critical temperature. 


3. Q Percent Aqueous GIvcol Solut i on Final Dpsipti - 

The critical aspect about the glycol coolant design was minimizing temperatures 
near the stagnation region. This was accomplished by using small tube sizes and as many 
tubes as allowed by the geometry. Since the density of this coolant is much larger than that 
of liquid hydrogen, and since the glycol flows in a closed system with a heat exchanger, 
high mass flow rates could also be used to minimize initial wall temperatures. This, 
however, created large pressure losses which would have to be overcome with a powerful 
pumping system. 

The following configuration was determined to best satisfy the design constraints 
on wall critical temperature and coolant pressure loss (see Appendix 3.2): 
mass flow rate : 150 kg/s 
Tube diameter : 6.5mm 
wall thickness : 2.5mm 
initial number of tubes : 20 
Branching locations : 0.1m, 0.75m, 1.9m, 4.35m 
The velocity of the coolant remains within the sonic limit both in the impingement 
nozzle (953 m/s) and the tubes (224 m/s). The outlet temperature of the glycol solution is 418 
K. The pressure loss in the system is relatively high (1422 psi). Therefore, all of the 
constraints were satisfied save that of pressure loss. 

This high pressure loss was found to be unavoidable. In order to keep the critical 
wall temperature below 1000 K at all stations, a high mass flow rate and a thick skin were 
required. Since weight is a critical issue in any aircraft design, the 2.5 mm wall 
thickness was judged to be as thick as desired. Taking all of the above into consideration a 
compromise on the pressure loss was required. This resulted in a pressure loss 722 psi 
greater than the imposed design constraint. 


Because the ethylene glycol design required*heat exchanger, a preliminary heat 
exchanger analysis was performed in which heat transfer coefficients of liquid hydrogen 
and glycol were compared. Given a tubular counter-flow heat exchanger, calculations 
were performed (listed in Appendix 1.5), and it was determined that the heat transfer 
coefficient for the liquid hydrogen was only about two to three times smaller than that of the 
ethylene glycol, given similar conditions. Thus, a configuration of bare concentric tubes 
could be used without fins on the inner tube extending out into the liquid hydrogen. If the 
proposed ethylene glycol design was chosen as the preferable design, the next step then 
would be to utilize software capable of designing such a heat exchanger. 


12 . Comparison of the Two Designs 


The liquid hydrogen design met all of the design constraints with a satisfactory 
margin of safety, whereas the ethylene glycol design failed to meet the pressure constraint, 
and met the temperature constraint with a very low margin of safety. In addition, the 
ethylene glycol design required a heat exchanger and a large pump, as well as a greater 
skin thickness, all of which add considerably to the weight of the drone vehicle. Therefore, 
the liquid hydrogen design was deemed the better of the two. 

One concern about the liquid hydrogen design is safety. In particular, a leak may 
prove to be disastrous for such a highly pressurized system of hydrogen. Because the will be 
placed on an unmanned drone vehicle, weight considerations dictated this small 
compromise in safety. In application to the NASP, however, appropriate safety measures 
could be taken to insure against such an fftturonoo. 


8 Summary 


Given an existing configuration, the authors determined a suitable flight condition for a 
hypersonic drone that would require some form of active cooling scheme and provide the 
propulsion group with an adequate challenge. The flight condition chosen was Mach 10 at 
100,000 feet. The authors wrote a FORTRAN program called HETAQ (Hypersonic Equilibrium 
Temperature and Heat Flux) that parametrically varied Mach and altitude while computing 
equilibrium temperatures, heat fluxes and flight times for each flight condition. 

Heating rates for regions aft of the stagnation point were derived for the chosen flight 
condition. An active cooling configuration was determined and consisted of coolant impinging 
on the inside of the skin at the stagnation region and then being routed through numerous tubes 
running back along the inside of the nose cone. 

The authors wrote ACSES (Actively Cooled Structure Evaluation System), a FORTRAN 
program that allows the user to vary tube and manifold geometry, mass flow rate, and coolant to 

observe what effect these parameters have on critical temperature, coolant temperature, and 
pressure loss. 

A design for liquid hydrogen coolant was formulated that optimized pressure loss and 
mass flow rate, given imposed design constraints. A design for ethylene glycol was also 
formulated. The hydrogen design was deemed the better of the two for the drone vehicle, based 
on weight considerations. 
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Skin Equilibrium 
Temperature Calculation 

Energy balance on skin assuming no 
heat transfer through back-face: 


convective heat _ radiative heat 
transfer ~ transfer 




ZCA 




q s = heat flux at stagnation point from 
Lees approximate method 

A = area 

£ = emissivity of material 

c = Boltzmann constant 

T v eq = equilibrium wall temperature 
T m — ambient air temperature 

T was calculated as a function of 

^ gq 

altitude, Mach number and nose 
radius to obtain a data base 



EQUILIBRIUM TEMP, (deg R) 


MACH NUMBER 

~00 8.00 9.00 10.00 11.00 12.00 13.00 



o 


Y- • 



Approximate energy balance to give an 
estimate of the flight time required to 
achieve the equilibrium skin temp. 

p CV(T v eq - ) = 

q s AAt-aGA(T veq .T.)At 

Some sample materials 
at 100,000 ft; Mach 9; radius =1" 


flight time =1.8 min 



Skin Thickness 


1 mm 
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Aluminum: 
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Titanium: 

2.7 sec 
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* Maximum skin temperature is midway between coolant tubes 


* Tubes welded to skin chosen to simplify fabrication 






illDinq ernent cooling at noss 



• Livingood's correlation for Nu = f(Re, geometry) 

Nu- 2,98 (Re) 0 585 f £ ,1.10 m -0.007 

v D ; 1 d ) 

D = hemisphere diameter 

zn = nozzle exit - to - target separation 

d = nozzle diameter 


• Average heat transfer coefficient. 


• (Tin - Tnoz) = 


• • 

q - q 

~ nose 2 rad 

area * h 


h = Nu * k 
d 


h>% Z 


Appendix 1: Sample calculations 



Appendix 1.1: Stagnation Point Heating 
(Lees' Approximate Method) 
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Appendix 1.2: Heating Rate Away From Stagnation Point 
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Appendix 1.3: Heating rates for Noses of Various Radii 
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Appendix 1.4: Equilibrium Temperature Calculation 
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Appendix 1.5: Preliminary Heat Exchanger Analysis 
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Appendix 2: Computer Codes 


Appendix 2.1.1: HETAQ Description 



HETAQ Manual: 


We used Lees Approximate Method in HETAQ because it quickly and accurately deter- 
mined stagnation point heating rates. From there, obtaining the skin equilibrium tempera- 

tures for varying flight conditions was just a matter of equating two equations and solving 
for the unknown. 

Using the simple energy balance equation: 
convective heat transfer = radiative heat transfer 
q * = q T 

the convective heat transfer is identical to the radiative heat transfer So we can 
now equate q, with the radiative heat expression: 

q, = q T = ea(T4 all - 7*) 

Now, equating the expression of ,, from Lees Approximate Method with that from the 

radiative expression and solving for the unknown wall temperature T wM gives a fourth order 

equation with which Newton's Method was used, iterating to solve T„„„ with the HETAQ 
code. 

Due to the limitations of fuel it was necessary to determine how long it would take for 
the nose to reach its equilibrium temperature. If it is longer than the calculated flight time, 
then the data obtained is in the transient region, and possibly not as valuable as data 


taken at a cruise equilibrium temperature. Therefore, this time calculation was a valuable 
step in defreffifti a g what the proposed flight condition would be. 

In calculating the time it took for the skin to reach this equilibrium temperature, an 
approximate energy balance was used to give an estimate of the flight time required: 

pCV(T w , e q — Tinitial) ~ q s AAt — G£A(T w , e q — T 00 )A,t 

This equation assumes that the nose piece is a simple thin skin structure (no heat sink) 
which is actively cooled to the desired wall temperature-to be determined by the material’s 
thermal properties. The time needed to reach equilibrium is quite short for the materials 
we analyzed. Using a thin titanium skin around 1 mm thick would only require* time of 
1.8 minutes. This should leave enough time for equilibrium data to be taken. 

Here is a summary of the process used. Calculate what the stagnation point heating 
rates will be using Lees Approximate Method. Then, determine what the skin equilibrium 
temperature will be, using the q s = q T energy balance equation. Now, in order to show that 
there will be a heating problem, if the equilibrium temperature is greater than the limiting 
temperature of the nose material, active cooling must be used. To determine how muchj 
cooling will be necessary, set the wall temperature at the limiting material temperature and 
use Lees Approximate Method to determine 

Qa — Qr Qcool 

where q coo i is the cooling rate to be determined by the active cooling system. 
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Appendix 2.1.3: HETAQ Results 
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•LIST 


1 ° SP ( 12 ) ' TS ( 4 ) 12 < 6 ' 4?’ 6),FL<4 ' l2,6) ' Ti?(4 ’ 12,6) ’ Ti<4 ’ l2 ’ 6),T2(4 ' 12,65,1 

2 Tp(T : HOME : READ N: FOR I = 1 TO N: READ Z< I ) ,RO< I ) ,TI( I) ,CS( I ) s NEXT 

I: RlM READ IN DATA... ALT (FT). DENSITY (SLUGS/FT‘3) , TEMP INFINITY 
(R). SPEED OF SOUND (FT/SEC) 

4 READ NUM: FOR I = 1 TO NUM: READ M$( I ) ,CP( I ) . DNS( I ) .THC I ) : NEXT : REM 
__ RtAD iN MATERIAL PROPERTIES 

5 FOR I = 8 TO 12: READ ISPd): NEXT : REM ISP FOR MACH # 

6 E ~ 3 - 141592654: S° = 605: REM SET SIGMA. EMMI 

SSiWTr, Pi , PlAN FORM ARlA (FT'2) 

7 Z = 0.01: REM ACCURACY SETTING 

8 DI = 8 : L = 40 :CD = 0.20:HD = 4.42: REM set dimensions of tank, draa co 

density of liq hydrogen 

9 V °<.. = -a PI r # (<DI / 2) 2) * L:W H = HD # VOL: REM voi , weiaht of H2 in 

r t 3, LO 

10 PRINT : PRINT ‘ ( A)utomat ic calculation or <M)anuaI' > INPUT T 

*: IF IS = "M" OR 1$ = “m" THEN 2000 

97 REM 

98 REM VARY PARAMETERS 

99 REM 

l"' FOR I = 4 TO 4 

1- . FOR M = 8 TO 12: REM MACH NUMBER 

• FOP K = 1 TO 6: REM RADIUS IN INCHES 

2- - R = SOP CK / 12) 

210 CS = CSC I) 

220 TI = TI( I ) 

230 RO = RO( I) 

240 ISP = ISPCM) 

500 GOSUB 5000: GOSUB i 0000 : GOSUB 30000 

= INT u0 ° * T) ' 100:0c I.M.K) = I NT c 100 * 0) / 100 

610 h Z 1 T ° NUK:7S( 1 *M,K.L> = INT (100 * DNS(L) * CPCL) * THCL) * 

CTCI.M.iO - 500) / Q(I.M.K)) / 100: NEXT L 

620 F n-u’ M Lp = (i0 ° * FL) 7 100 :TRC I ,K .K> = INT <100 * TR) / 100: 

_RaM ul= FUEL FLOW. TR = THRUST 
bd 7: C I.M.K) = Ti:T2( I.M.K) = T2 


997 

REM 


998 

REM 

SAVE D 

999 

REM 


1000 

PRINT 

CHRS 

1010 

PRINT 

CHR$ 

1050 

PRINT 

" I n 1 1 i 

11 

: HD : ” 

lb/ft ‘3 

= 

" :DI ; 1 

' ft" 

1060 

PRINT 

"Fuel : 

1070 

PRINT 

"Cd = 

: ] 

PRINT 


1 1 "O 

FOR I 

= 4 TO 


- - . * i. uj ubinu . un i n 

1050 PRINT "Initial conditions:": PRINT : PRINT "Density of hydrogen = 

Spec 3): Tan*: length = 11 :L:“ ft Diameter 


area = " :S0:" ft '2 


PRINT : PRINT : PRINT “Altitude = " : 2 < I ) : “ ft' 
FOP M = 8 TO 12 


1200 
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_ 1 ') PRINT SPCc 3) j ■' Mach = : M 

louO FOP K = 1 TO 6 

] PRINT SPC( 6): "Radius = : ‘;Kr inches" 
i <00 PRINT SPCc 9): "Nose egi»i I temp = " i T ( I , M , K 3 1 v R' i Spec :i)rq = « j 
Q(I.M.K);" otu / ft '2 sec 

1710 PRINT : PRINT SPCc 9):"Steaav state reached 1 mm 1/ 

16" : CHR$ (34);" (1.59 mm)" 

— 1720 PRINT SPC( 20) ;MSC 1 ) : “ : SPCC 5> :TS< I ,M,K, 1 > ;” sec"; SPCc 8);TS 

( I , M , K , 2 ) ; 11 sec" 

1730 PRINT SPC( 19) ;M$( 3) : " : SPCC 5) ;TS( I ,M,K,3> ; “ sec"; SPCC 8);TS 

(I.M.K.4):" sec" 

“■ 1800 REM FLIGHT CONDITIONS 

1810 PRINT : PRINT SPCc 9);*Fiight conditions: ISP = “ ; ISP(M) ; " sec 

Thrust = 11 ;TR( I ,M , K ) : “ Ids 1 ' 

— 1820 PRINT SPCC 29) : *• Fuel flow rate = “ :FL( I ,M,K) : " io/sec": PRINT SPCC 

29) ;“F1 lght time = " : T 1 ( I , M , K > : “ sec Flight time = " :T2C I ,M,K) : “ m 
in" 

___ 1890 NEXT K.M.I 

1900 PRINT CHRs (4): "CLOSE HETAQ. DATA 

1901 END 
1997 REK 

— 1998 REM DO MANUAL CALCULATIONS 

2000 REM 

2010 PRINT : PRINT : INPUT "Mach numoer = " :M 
___ 2050 INPUT "Density = " :RO 

?~"0 INPUT “Speed of sounc = :CS 

2 0 I NPUT l! T( infinity) = 1 ; T I 

?' '0 INPUT “Nose radius = (in ft)“ :K:R = 3QR ( K ) 

— 2 0 ISP = ISP(M) 

2200 GOSUB 5000: G0SU5 10000: GOSUB 30000 

2205 PRINT "Initial conditions:": PRINT : PRINT "Density of nvaroaen = 

— J : HD : " lD/ft'3": SPCc 3):"TanK: length = ":L:“ ft Diameter = " : D I ; “ 
f t " 

2206 PRINT "Fuel volume = ":V0L:“ ft *3"; SPCC 5):“Fuei weight = " :WH 

2207 PRINT "Ca = ‘ : CD : SPCc 5>:"Pian form area = ";S0:" ft‘2": PRINT 

2210 FOR L = 1 TO NUMiTM(L) = INT c 100 # DNS(L) * CP(L) * Th'(L) * CT - 

500) / 0) / 100: NEXT L 

2300 PRINT : PRINT : PRINT * Equ i i i or l um temp = “:T:“ R“ : PRINT “Heatina 

— rate = " :0:‘ Dtu / ft'2 sec 

' : 9>: Steady state reached i mm 1/ 

16": CHR$ C 34 ) : 1 <1.59 mm)' 

2320 PRINT SPCC 2C):MSC1):": ; SPCC 5):TM(1):" sec"; SPCC 8):TM(2):" 

sec" 

2330 PRINT SPCC 19 ) ; MSC 3) ; " : " : SPCC 5 ) : TMC 3 ) ; sec": SPC( 8):TMC4);“ 

sec" 

— 2400 REM FLIGHT CONDITIONS 

2410 PRINT : PRINT SPCC 9); "Flight conditions: ISP = ":ISP:" sec Th 

rust = ";TR;“ Ids" 

_ 2420 PRINT SPC( 29):"Fuei flow rate = " ; FL : " lb/sec“: PRINT SPCC 29); 

“Flight time = “;Tl;“ sec F i i ah t time = ":T2:" min" 

2500 GOTO 2010 
- 7 REM 

— 4990 REM CALCULATE COEFFICIENTS 

‘ * REM 

5000 V = M * CS:A = (2.192031E - 8) * ( SOP C RO ) ) # (V * V * V) / <E * S 
__ ) 

5010 B = C T I / 4) ♦ CV # V) / 50100 :C = (A / B) * IE - 5 
o020 D = (A / B) * 0.234 

5030 REM M = MACH #: CS = SPEED OF SOUND ( FT/SEC): V = VINFINITY (FT/S 
EC): RO = DENS I i x AT ALTITUDE: TI = TEMP INFINITY (R): S [=] BTU/FT* 
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2 SEC R' 4 
5040 RETURN 

9998 REM 

9999 REM DO ITERATION CALCULATION 

10000 REM 
10100 T = 2500 

10500 X=T- <<R *T * T #7 * T) + (C # T # T) + (D * T) -A) / ((4 # 
R*7*T*T)+(2*C*T#7)+D> 

10800 IF ( ABS <T - X) > Z) THEN T = X: GOTO 10500 
10900 Q=E*S*T#T*T*T 
11000 RETURN 

29999 REM 

30000 REM DO THRUST CALCULATIONS 

30001 REM 


30010 TR = CD # 0.5 * RO # V # V * SO: REM THRUST EON 

30020 FL = TR / ISP: REM FLOW RATE (FT '3 / SEC) 

30050 Tl = INT (100 * WH / FL> / 100 : T2 = INT (100 * Tl / 60) / 100: REM 

FLIGHT TIME IN SEC AND MIN 


31000 

RETURN 

49997 

REM 


49998 

REM 

ATMOSPHERIC DATA 

J-OQ9 

i .00 

REM 


DATA 

5.70000,1 .3761E-4, 392. 97. 971 

r ' ' 10 

DATA 

80000.8. 445E-5 , 398 . 46,977. 

5 20 

DATA 

90000,5 . 21 73E-5 , 403 . 95 , 984 

50030 

DATA 

1 00000, 3. 320E-5, 409. 17,991 . 

50040 

DATA 

110000. 2. OlE-5. 419. 87. 1003.; 

60000 

DATA 

4, "Aluminum 1 ,0.208.169,3.28 

60010 

-3 

61000 

DATA 

"Magnesium" .0.232.109.3.281; 

DATA 

2350.1988.1775.1488. 13i3 


,0 .208 , 169,5 . 208E 
,0.232, 109.5.208E 
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This program integrates the theories of conduction, convection, and fin 
analysis discussed already. Its main objective is to provide a first approximation 
in the design of a combined impingement and actively cooled panel concept for the 
drone nose. The code allows the user to vary tube /nose geometry and 
material/coolant properties with which it uses to calculate the critical skin 
temperature distribution, total pressure loss, and exiting coolant temperature. 

The design constraint is that the skin of the nose never exceed a maximum 
allowable temperature. A branching routine is incorporated in the program to 
nelo avoid this, but at the expense of an appreciable pressure loss. Thus the 
task is to find a design that limits the abount of branching necessary to keep the 
skin temperature deiow this maximum value and the pressure losses at a minimum. 

ACSES reauires the following input parameters: 

- Numoer of cooiant tubes (tubes) 

- Maximum distance back on nose to be analyzed [m] 

- Number of stations (subdivisions) for iteration of properties (xnum) 

- Coolant tube inner diameter (dtube [mm]) 

- Thickness of coolant tubes (xtube [mm]) 

- Friction factor of coolant tubes <f> 

- Radius of the nose cone (rnose [mm]) 

- Skin thickness of nose cone (xskin [mm]) 

- Skin limiting temperature (tmax [O 

- Mass flow rate of coolant into the manifold (mfrO [kg/s]) 

- Impingement nozzle inner diameter (nozd [mm]) 

- Distance from nozzie exit to impingement on back of nose wall (zn [mm]) 



- CJhMitfs Hi (u«mn HvHr-HHBfi, water invent sHiutiwfi caB**), ThermihHt $0) 

- Initial coolant temperature (tman [O 

ACSES gives the following as output: 

- The location of each required brancn 

- At each station along the nose cone: 

— Surface heat flux (qmm [W/m‘2D> 

— Bulk mixing temperature of coolant (tmav Ck]) 

— Inside tube wall temperature (tw CO 

— Base of tube temperature (tO CK]) 

— Local Revnolds number 

— Pressure loss (pressl CPaD> 

- Final results 

— Exiting coolant temperature (tmo CO 

— Total pressure loss (ploss [PAD, [psiD, Catm3> 

— Outlet tube mass flow rate (mfr Ckg/s3) 

— Total number of branches 

Assumptions: 

The properties of tne tube and skin materials were evaluated near their 
maximum allowable temperatures, tmax. and taken as constant along the entire 
nose cone region. Since the properties don't vary too greatly at these 
temperatures, a constant value serves as a good approximation. 

- Nonaccelerating, level-cruise at Mach 10, 100,000 feet. This flight profile 

allows ACSES to calculate the stagnation point heating rate and use it to 
determine the heat flux at any point along the nose. 
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- Each tune branch is perfect and occurs over a negligible distance. Pressure and 

frictional losses are neglected as well. 

- Coolant properties vary greatlv along entire nose region but are taken to be 

constant over each differential station. They are evaluated at the mean bulk 
temperature (Tmo - Tmi) / 2 and iterated to within 1 degree accuracy. 

- Tubes start at hemisphere of nose cap and extend back radially. 

- The distance that the heat transfer across the coolant tubes occurs is taken to 

be 50% of the tube circumference. This approximation is valid due to the fin 
effect of the tube shell. Analysis shows that the wall temperature decreases 
to the coolant temperature within 15% of the circumference. The other 35% of 
the circumference approximates the heat transfer through the solder that 
welds the tube to the nose skin. 

- Radiation off me nose and awav from the inside wall have been neglected. 

Calculations show mat it contributes less than 1% of the total heat flux. 

Eauations: 

After the inout of the parameters, the program goes through a set of 
eauations to calculate the output. The order is as follows: 

- Find the coolant temperature exiting the manifold 

- Calculate nose and wall enthalpy to be used in heat flux calculation 

h = 

- Start at first station 
— Calculate heat flux 

q = 
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— Calculate distance between two tubes 
L = 

Iterate on energy balance to find temperature leaving station 

Evaluate Cd at Tm 

— Find temoerature leaving station 
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Tmo = 


Recalculate mean oulk temperatue 
Tm = Tm] t Tmo / 2 

— Check with first guess, if not within 1 degree, iterate again 
Evaluate coolant properties at new mean bulk temperature 
Calculate Prandtl. Reynolds, and Nusselt numbers; 

Pr = 

Re = 

Nu = 

from turbulent inner flow correlation valid for Re greater than and 0.7 < Pr < 

13000 ano Visrositv is evaiuateo at mean bulk temperature and at wall 

temperature. 


-- Calculate heat transfer coefficient between inner wall and coolant 


h = 


— Calculate thermal resistance network 


R 1 = 


R 2 = 
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— Wall temperature of tube from convection of coolant is used to evaluate 
viscosity of coolant — iterate to within 1 degree as before. 

tw * tmav (am * Rl) 

— Find base of fin temperature from conduction through tube wall 
tO = tw + (qm * R2 j 

— Find critical wall temperature using the base of fin temperature with the fin 
analvsis 

Tcrit = tO * (amm * (L ' 2) / (8 * kskin * xskin)) 

— If tcrit > tmax branch 

— douole the numoer of tubes 

— half the mass flow rate 

— Evaluate pressure loss along the station. Total pressure loss is found from 
the sum of each station loss. 

ploss = 

— Print out calculations to an output file 

— Continue the analvsis with the next station 
- Print out results 

— Total pressure loss 

— Final coolant exiting temperature 



Limitations of ACSES: 
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The program only permits the design of a welded tube system. This 
configuration was chosen because of its ease to manufacture. The analysis also 
assumes a uniform heat distribution along the skin which In reality, might have 
several "hot soots" due to surface irregularities. This could be potentially 
dlsasterous because the existing design has no way of knowing where these hot 
soots are or anv wav to divert coolant to cool them off. If the temperature 
exceeded critical limits, the skin loses its structural integretv. This in turn 
would cause the cooling svstem to fail and allow temperatures to reach 
eauilibnum and burn through the skin, destroying the drone. Currently under 
investigation are wavs of discouraging tne formation of hot spots — possiblv bv 
aoolving a tmn ceramic layer on the outside of the nose. 

Advantages of ACSES: 

The orogram is a Quick and relatively simple wav to initiate the design of an 
actively cooled nose cone. The assumptions and correlations used in the code are 
accurate for me Purposes they are sprvmg and were easily integrated into the 
analysis. ACSES Drovioes the means for a first approximation on what might 
otherwise be a series of tedious calculations. Due to the interactive nature of 
the program, many different cases varying in geometry, coolant choice, and 
materials can be run and compared in a short amount of time. 

In the future. ACSES can be modified to provide the thermal design of the 
leading edge of the drone as well. Correlations would be needed for the given 
wing geometry which exDress the heat flux along the leading edge. The code 
would also have to allow for variable mass flow rates with multiple 
implndgements along the span. 



Appendix 2.2.2: ACSES Listing 
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fTfldl Code 

s 



c 


A C S E S s 
Ginger Koch 


Actively Cooled Structure Evaluation System 
James Marco lesco , Greg Noffz? Jerry Yen 


c 

c 
c ■ 
c 
c 
c 


Vers i on : 5 . 2 A . bS 


INPUT PARAMETERS : 


c 

tubes 

c 

d tube 

c 

x t la b e 

c 

x s k x n 

c 

x num 

c 

inf r 0 

c 

no zd 

c 

choice 

c 

p r n t 

c 



u. 


- number of tubes 

= tube inner diameter? input Cimn] , converted to Cm] 

— t h i c k n ess of t u be? i n p u t C m m3? co n v e r ted to C m ] 

= skin thickness? input Emm!]? converted to Cm3 

= number of stations 

-- rn ass flow rate of c o o 1 an t i n to m a n i f o 1 d C k Q / s 3 
i mp i n a e m e n t nozzle i n n e r d i a m e t e r C m 3 
= c h o i c e o f c o o 1 a n t 
=-- 1 ) shor t pr i nt 2 ) 1 ong pr i r« t 


c 

c 

c 

I 

L. 

c 

c 


c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


C 0 N S T A N T F ARAMETERS s 


r n o s e 
a .1 p h a 
h 

z n 

tmax 


tmaxr 
t rn i a n 

Pi 

f 

ksk i n 
k tube 


-= r a d x u s o f no s e c o ne C m 3 

•--- cone angle? convert from degrees to radians 
= m a ' i m u i n cl i s t a n c e b a c 1: f r o m n o s e a n a 1 y z e d [ m 3 
— i m p i nq em e n t no z z 1 e d x s t a n c: e E rn 1 
=" mater iai 1 xml ting temper ature EK3 
1000 Ck3 for molybdenum 
~ IS 0 0 L R II f o r m o 1 y b d e n u m 

-• t emp er atu r e o f c o o lant i n t o m a n i t o 1 d C h 3 
“■ 3.1 9 1 5 9 2 & 5 9 

= f r i c t i o n fact o r of c o o 1 a n t t u b e s 

t hi e r m a 1 c o n <3 u c t i v i t y c* f n o s e C W / m K II 
e v aluatsd at t m a x f o r m a t e r i a 1 
i •: j 3.77 f o r m olybde n u m at T m a x 
thermal conductivity of the tube wall EW/mK3 
= 2 9 . 4 [ W / m K 3 f o r C u p r o n i c k e 1 3 0 % 


Comrner c i a ] Tub i nq S l zes s 


c 

ID 

GD 

Th i c knees 

gage 

max nui 

c 

[mm3 

[ i n 3 

[mm3 


for 1 " 

c 

c 

1 .63 

1/3 

0 . 76 

A 

50 

c 

3.29 

3/16 

0 . 76 

A 

33 

r 

9.83 

1/9 

0 , 7 6 

A 

25 

f 

3.92 

1 / 9 

1 . 29 

CO 

25 

1 

6 . 32 

5/16 

0 . 8 1 

2 1 A 

20 


c 

c 

c 

c: 

c 

c 


Cupronickel 70% Cu 3C*/« Nx Wrought 
y i e 1 d s t r eng t h : 57 ? 000 C p s i 3 393 C MP a 3 
tensile s t e nq t h : 60 ? 00 0 C p s i 3 919 C MP a 3 
m e 1 1 i r i q t e mi p e r a t u r e : 15 0 0 C K 3 

so ec; l f i q r a v x t y r. 8.95 


Tubing Sizes: 


OD Thickness cjaqe 

C i n 1 II mm 1 


1/2 0.81 21 A 

1/2 1.65 16 

5/8 0.89 20 A 

3/ A 0.89 20 A 
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N D E N T P ARAMETERS : 


= coolant temperature exiting manifold 
and entering tubes = tkl CK1 
rr t e mo or ary t f i 1 m cal c: u 1 a t i o n E K 3 
•■= bulk mixing temperature of coolant, used 
to evaluate properties CK3 
base of fin temperature EK3 

- temperature inside tube wall CK3 
~ m i d 1 i ne critical t emp er atu r e E K 1 

= m e a n f 1 u i d t e mi perat u re at static n k ? i n a n d o u t E K 3 

--- f i n a 1 fluid e x 1 1 i n q tern p e r a t u r e C K l 

- m a s s f 1 o w rate i. n o ne t u b e C k q/s] 

-- no z z 1 e velcci t y E rn / s 3 

= t u b e v e 1 o c i t y C m / s 3 

™ t o t a 1. h e a t at d x s u r f i n t e r m s o f e n t hi alp y f r o mi free 

s t v e a m c o nd 1 1 i o n s C W 3 
= CW/ml 

- [W/m*'*E3 

= c hi o i c e of c o o 1 a n t 

-■ density of coolant from curve fit Ckg/m "311 

- d y n ami i c vise o sit y o f c cola n t f r o mi c u r v e fit E kg / m v 3 .1 

” d v r i a m i c vise o sit y o f c o 1 1 a n t at wall t e rr? p 

= s p ec i f i c h e a t o f c o o 1 a n t f r o m c u r v e fit E J / k q K 1 

= p r a. r id 1 1 n u m b e r o f c o o la n t f r o m c u r v e fit 
= : = P e / n o Ids n u m b e r of c o o lan t 

“ N u s s e 3. t n u mbe r o f c o o 1 a n t f r o m c o r r e 1 a t i o n 
= n o s e e n t h a 1 p y 

= w a 1. 1 e n t h a 1 p y 

= local heat transfer coefficient for turbulent flow 
c o r r e 1 a t i o n at static n x C W / m ' ' £ K 3 
= the r m a 1 c o n duct l v 1 1 y o f f 1 u id f r o rn c: u r v e fit C W / mK 1 

- system pressure loss C?J 

= dummy variable? piecewise pressure loss 
= distance between tubes at distance x Cm3 
= position number? k = E? xnum 

= distance back from manifold of next and previous 
station = < k — 1 ) ■* h / x nu m C mi 3 
= distance between stations < xE - xl) Cm3 
= area that heat transfer takes place across tube 
= area that conduction takes place across skin 

= number of branches 
= distance back of each branch 
= flag that branching has occurred 





( 

L. 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c: 

c- 

c 

c 

c 

c 
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COOLANT PROPERTY SUBROUTINES: 

thcond (choice, tmav.kflu id) thermal cond at tfilm, kh CW/mKD 

specht (choice, tmav,cp ) specific heat at tfilm, cp C J/kgK3 

dynv ( choice » tinav < inu ' dynamic viscosity at tfilm, tw 

mu E k g / ( m ) < sec ) 3 

d e nse (choice? t mav ?ro) d e ns i t y ^ t tfilm? e x p o nen 1 1 a 1 fit? 

ro E kg / m % 3 II 

fit ( c h o i c e ? t ro a v ? a ? a n s w er ) c u r v e fit r out i n e s 


M A I N P R G G R A li s 


c 

C234567 


( 


p r o q r am m a i n 

1 n t e q e r x 1 1 u m ? tub e t? f 1 a q ? n u m tub ? c h o i c e ? n b r ? b r 
real hi mt ro ? nrfr ■. kfluxi d ? k sk i n ? k tube ? h x ? nu «< L * hnose < liwa 1 1 ? f ? 
no zd ? nve ]. ? d t ube * mu ? mow ? xhr ? rna x 


c 

c R E A D I N 1 N R U T D A T A : 

c 




open ( 2 ? f i 1 e - ? ac:ses . i n ^ ) 

r e a d ( E ? 1-4 O 0 ) r n ose?alpha? z n 

1 400 f c« r mi a t (4 f 1 0 - 5 ) 

i ead , E i 1 40 1 ) t met x ? t met x r ? k s J- ; i n ? ero l ss ? k t ube ? ep) s 1 o n 

1 40 1 forma t ( 6f 1 0 - 3 ) 
e 1 ose ( E ) 

pi = 3.141592654 
pi oss = 0. 

1 00 format ( / ) 

200 format ( // ) 

print * , Parameter input:* 
wr i te < 0 ? 1 00 ) 

print * ? * Coo la n ts: L 1 J liquid h y d r o q e n ? 

print * ? * E23 water-glycol solution ( 30 * 4 ) " 

p r i n t * ? E 3 ] Then m i no 1 60 * 

print *■ ? Input coolant choice > ■ 

read *? choice 

P v " i n t * *» Mi n i m u m c o o 1 a n t t e rnp e r aturess ’ 
pr int *?•■’ E I 3 12 K 

print * ? - EE 3 270 K 

print * ? * E 3 3 250 K 

P r i n t ■*- •» I np li t 1 n i t i a 1 c o o 1 a n t te mp E f : ' 3 > =• 


IOC 


re a a 



h t ft : 3. i “( 

print * , • Input tube diameter C mm 1 

read *• n d tube 

cl t ube = dtube / 1000 

print * i Input tube thickness [mm3 

read * % x t u b e 
x t u b e = x t u be / 1 C> 0 0 

max = 2 * pi * rnose / (dtube + <E*xtube) ) 
print * 9 For given geometry * # tubes < J jiriax 

print * ? ’ Input number of tubes 

read *• * tubes 



p r i n t * ■« * I npu t wa 1 1 t h i c k ness C mm 3 

read * * xsk i n 

:: Si: 3. i i =■ X £ k i H / 1000 


print * n ’ Mo.;- distance back on cone Cm3 



h nose — ( 0 . 83m * tmarr ) +• (0.00001 * t max r ** Li. > 
hwall - hit use 
r no s e i = r no s e * 3.2 308 

c: 1 -"-9 792* ( r nose :L * *■ ( - . 3 > ) *a 1 pha * ( 2090 . -hnose ) 

* / k 8 c* tv ■ . ■“* > 3 . 6ot V *l i w a 1 1 ) 



F I N D C 0 0 L A N T TEN P 
E X I T I N G I M P I N G E hi E N T : 


c 


/<?/ 


e d s 1 on ) 


c 


tiTu = tcin 
tmav = 1.2 * tm i 
tw = 1 0 . * t m i 

l - o. ORIGINAL PAGE IS 

d x = h / xriUiTi 0F POOR QUALITY 

b r ~ 0 
nbr = 0 

num tub = tubes 


i 

a* 

a 

= 

alpha 

*• 

E * | 

pi / 360 


f - 

0 

► . oo a 





mfr 

= 

mf rO 

/ 

num 1 1 

Lib 


xh - 


pi * 

d t 

iibe / 

C. M 


XS - 


* 

Q. 

d 1 1 

u b e / 

ci « 


cal 1 


dense ( 

<:: h 

o i c e , 

t ffi i i r o ) 


nve 1 


=- (A*fr 

f r 1 

0 > / 

y r o * p i * 

(nc« :-?d ** E)> 

t v e 1 


-- ( A* ill 

if r 

) / 

( r c* * p i * 

(dtube ** 2) 

o p evi 

( 

5 , f i 1 e 

=s 

* acs< 

as • out * ) 


o p e n 


1 0 , f i 1 

e 

“ 'pi 

o t . q s ) 


o p en 

( 

1 1 , f i 1 


- ' p i < 

:« t . tut > 


o p e n 

i \ 

1 2 , f i 1 

e 

» ’pi 

o t . t c v ) 


wr 1 1 

e 

( 5 , 1 0 

('8 




wr i t 

e 

( 5 , 1 0 

09 

> 



if ( 

o 

v n L „ eq 

. 1 

) thei 

n 



wr i te < * i 1 008 / 


Wi" 1 1 1-.‘ n j. OO'r ) 

end i f 


^ S T A T I Q N C A L C U L. A 7 I 0 N S : 

c 

c iti a. i n rout i ne to c ci 1 c u 1 8 1 1 ? troi <. tmo at esch (ox) s t s t i u t 

c to a. d i stance h away from manifold’. x r.u m times., 

c ]. i in i t s o t i n t e o r a t i o n a r e x 1 , x 2 

c 


do 10 k - ct n x num 
X 2 z - h * ox 


c 


c 

c 


c 


1 0 


L == (E. * pi / num tub ) * ( x2 * sin (rad) + rnose) 

C ci n v e r s i o n f a. c t o r f r o m B t j / s f t s q r d t o W a 1 t s / m s q r d 
is 11 358.8855 


qmm = 11352.8855* cl * (;■:£ *■* (-0,2)) 
qm = L * qmrn 
q = qm * dx 

if ( p r n t . eq . *0 ) wr i be ( * .. 1 5 1 0 ) k •. x 2 
for ma tOXO I3o’ x = % F 5 . 2 ) 

if (hr.eq.l) then 
b r == 0 

l f ( p I n t .eq - 0 ) wr i te (*,1011) nb r 
end 3 f 



c 


ITER A T I □ N 



c 


600 


700 


call specht (c ho i c e « t mav < c p ) 


tftio = tmi + <q / (mfr * cp > ) 


titer •- ( tmo ■+■ tmi) / 8 

i f ( ab=. ( t i ter - tmav ) . g t . eps 1 on ) then 
tuiav ~ titer 
go to 600 
end i f 


c a 11 t !' i c o nd < c h o i c: e * t nia v , k f 1 u i d ) 


ca 

I 1 

d y 

nv ■: ch 

i O 1 C € 

•. t m a v 

* mu ) 


c a 

1 1 

dv 

nv ( ch 

o i c: e 

* tw . m 

u w ) 


ca 

1 1 

de 

nse ( c 

ho i c 

e n tma 

v n ro 

) 

pr 

— 

mu 

* cp 

/ k 

f 1 u i d 



re 

- 

4 

* ffifr 

/ ( 

pi * 

mu #• 

d t u b e ) 

nu 


\ 0 

.. 087* 

e*-* 

0 . 8 > * 

( ( iTft.t .■ 

•' mu w )■*■*. 1 4 ) •* ( p r 

h x 


kf 

1m i d 

•* i»U 

/ dt 

ube 



r 1 = 1 

/ 

( b * 

h > 




1 1 j 1 1 e r 

“ 

tmav + 

( qrfi 

* 

r 1 > 


if (ab 

S < t; 

wx ter 

— t w ) „ 

c t; 

. ep s ]. c 

< n ) 1 1 1 e n 

tW - 

t V 

i t e r 





G O t O 

*7 






end :i. f 







r 8 = r 

t; j b 

e / i k 

t u b e •* 


h ) 


CO = 

t w 

+■ ( qm 

* r 2 ) 




t c r i t 

= t; 

0 -- { q i 

iilli * ( 

L 

•*** 8 / 

/ ( S * k s kin * xsk i n 

cress 1 


( 8*f *d 

* < m f r 

«*•* 

■a ) ) ' ■: 

< p i *• * 8 ) * r o * ( d t u be*-* 5 ) 

p 1 OSS 

=r P 

I o s s + 

press 

1 




•i f ( ter 1 1 . qe . < 0 . 9E * tma x ) ) then 
br = 1 

nbr = nbr + 1 
xbr = x8 

numtub = tubes * <8 *--* nbr) 
ffifr = (nfrO / numtub 
e nd i f 


1 00 b 

*■ 


1 009 


* 


format: ( • # 

’ tw: 

format ( " 

•- CK3 

w't :l te 1010 ) 


x : qs: t m a v : p r ess: 

t0 : tent: Re : ' ) 

Cm3 CW/mr 23 C K 3 CPa3 

CK3 CK3 ’ ,/> 

k n x 2 i qmm ? t ma v „ p r ess 1 jtw^t 0 ? ter it rre 


(03 


1 0 1 5 


1 0 1 o 

•* 


Wl 1 ^ < 1 o 9 i U 1 Zj > V Ct H Q (ft (h 

wr ite (11,1015) >:8 , tmav 
write (18,1015) >;8, ter it 
format (ei5.4,lx,el5.4) 

if (prnt.eq. 1) then 

write (*,1010) k , >;8 , qmm , tmav , press 1 
end i f 




j tWntOjitcriti. re 


format ( , 1 x , f 5 . 3 , 1 x , e8 . 3 , 1 * , eS . 3 , 1 x , eS . 3 , 1 x , 

eS . 3 , 1 x , e8 . 3 , 1 , eS . 3 , 1 , e8 . 3 > 

if ( br . eq . 1 ) th en 
write (5,1011) nbr 


101 1 


if (prnt.eq.l) write (*,1011) nb r 

format ( / , > branch # ^i^,*'’ < y 9 / ) 


end i f 

tmi == t mo 

v. 1 = ci 


1 0 continue 


c 

c 


F R I N T OUT RESULTS: 



wr i te < 5 .. 1 00 ) 
wr ite ( 5 , 1 090 ) 


•if ( prnt - eq k 1 ) then 
wr ite ( * , 1 00 ) 
write (*,1080> 
end i f 

1 020 f o r ma t { / / / / / / / / / / / / / / / / / / / , 

i f ( c h o ice. eq M 1 ) t h en 
wr ite (5,1 050 ) 


Resu I ts 


'"'d i n 1 1 1 a 1 pa r a me t er s 


1050 


1 060 


1 070 


1 f (ornt.eq.l) write (*,1050) 

format ( ’ Coolant considered s Liquid hvdroqen’ , /> 
el seif (choice.eq.8) then 
wr ite (5,1 060 ) 

if ( p r n t . eq „ .1 ) write ( * , 1 060 ) 

format < ’ Coolant considered: Water /qlycol solutic 

else 

write (5,1070) 

if' ( p r n t . eq . 1 ) wr ite ( * , 1 070 ) 

format ( ’ Coolant considered: Therminol 60’,/) 

end if 


write (5,1110) tcin 
f o r mat ( : ' T emper a t u r e e >: i t i ng 
write (5,1180) tmo 


manifold = 


1 1 1 0 


, e 1 0 . 4 , ’ C K 3 •’ ) 




i 1 f - *• 1 

t c* r ma n 


- in a i -e ; x r j n y r, e m \ j e r a r, li r e 

— i 0 1 U ■ 4 j 

" L R J ; 


1130 

wr ire 

( 5 * 1 1 30 ) 

p 1 C* S S n ( p 1 O 5 S * 1 . 4 5 0 40- 

4), (p loss/ 101350) 

t or mat 

( ? 

System pressure loss 

— , fc? 1 0 • 4 ? 

’ C Pa D 


* 


e 1 0 . 4 n 1 

’ Cpsi 3 ’ , fS.S, ’ atm-’ > 



/ 


wr i te 

(5,1 140 ) 

tubes 



( 

0/ 

1 1 40 

f or mat 

< 

Initial number of tubes 

= •’ , i 4 ) 



w r i t- e 

(5,11 50 ) 

(dtube * 1000) 




1 50 

t o r rn a t 

( 3 

Tube diameter 

= 3 , f 5 . 2 , 3 

[mm] 3 ) 



write 

<5,11 60 ) 

( >( tube * 1000) 




1 160 

forma t 

( ■’ 

Tube th ickness 

= 3 ,f 5.2, 7 

C mm I 3 ) 



wr i te 

(5, 1170) 

( r nose * 1000) 




1 1 70 

format 

< - 

Nos© radius 

= 3 , f 4 . 1 n 51 

C mm 3 3 ) 



wr i te 

<5,11 SO ) 

(xskin * 1000) 




1180 

format 

( ’ 

W all t h i c k ne s s 

= 4f4.1r 

C men 3 3 ) 



wr i te 

(5,11 90 > 

alpha 




1190 

f o r rna t 

( =* 

C o n e half angle 

= 3 , f4. 1 , 3 

C deg 3 3 

) 


w r i t e 

(5,1 200 ) 

mf r 0 



1200 

f o r it i a t 

( 3 

N o zzle rn ass f 1 o w r a t e 

= 3 , f 1 9 . 4 , 3 

Ckg/s3 

3 ) 


w r i t e 

( 5 , 1 2o2 ) 

nve 1 



1202 

for ma t 

( ■ Ve 1 dc i t y o f c o o 1 an t i n no c z 1 e 

= 3 ,f8.3, 3 

[ m/s 3 3 

) 


wr i te 

<5,1 2 04 ) 

mf r 




1204 

forma t 

< Outlet tube mass flow rate 

= 3 f 7 . 4 , 3 

C kg/s 3 

3 ) 




1206 

1208 

1210 

1215 

1220 


w r i t e < 5 .. 1 206 ) t ve 1 

format < 7 Outlet vel of coolant in tubes 
wr i t e ( 5 , 1 2 08 ) ( z n 1000 ) 

f o r mat ( I m p i n g e me n t t a r q e t d i 1 ci n c e = 

wr i t e ( u i 181 o ) ( no zd * 1 000 ) 

format; < 3 Impingement nozzle diameter = 
write ('5 * .1 2 .1 5 ) nb r 
for m a t ( ■ N u rn b e r o f b r a n c h e s = 

w r i t e < 5 , 128 C> > t m a >; 

f o r m a t ( 3 L. i nu t i n q m a t e r i a 1 t e m p e r a t u r e = 
c 1 o s e ( 5 > 

L. 1 C 1 S © ( 1 7 

c 1 use (11) 
c 1 o s e ( 1 2 ) 


’ if8.3, 3 [ rn / s 3 ’ ) 

, f 4 . 1 , 3 C mm 3 3 ) 


, f 5 . 2 ■ 


C mm ] 3 ) 


] 

?f5.0, 3 UK 3 3 ) 


f ( p r n 

t.e-q.l) tl 

i © Ti 

wr i te 

( * , 1 1 j. 0 ) 

tc in 

wr i te 

<*, 1 1 so > 

t rn o 

w r i te 

(*,11 30 ) 

p 1 o s S i ( p J c« s s 

wr l te 

(*,11 40 > 

tubes 

wr i te 

(*,11 50 > 

(dtube * 1000) 

w r i t e 

(*-11 60 ) 

< r tube * 1000) 

w r i t e 

< * , 1 1 70 > 

( r no 58 •t- 1 00O ) 

wr i te 

(*,11 80 ) 

( xskin * 1000) 

wr i te 

( * , 1 :i. 90 i 

a 1 p h a 

w r it© 

( * , 1 BOO > 

mf r 0 

wr i te 

( * , 1 EOS ) 

n v © 1 

wr i te 

( * , 1 E04 ) 

mf r 

wr i t e 

( * , 1 E06 ) 

t vel 

w r i t e 

( * , 1 SoQ > 

(zn *■ 1000 ) 

wr i te 

(*, 1210) 

(nozd * 1000) 

wr i t e 

<*, 1215) 

nb r 

wr i te 

<*, 1220) 

t max 


* 1 . *450 e-4 ) , (pi oss / 1 0 1 350 ) 


end i f 
end 


C 


c 

c 

c 

c 

c 

c 


SUBROUTINE NANI 


Jerry Yen 

Variable descr iptio n 

Pi i sigma ... c o ns t a n t s 
theta . . . angle from no st 
manifold ( dea '> 


version 5.26.88 


e tip of semihemisphere to back end of 


\ot> 


L 

c 

c 

c 

c 

r-' 

u 

. r 

r 

C 

C 

c 

c 

c 

c 

c 

c 

c 

c 

’ c 
c 
c 
c 
c 
c 
c 
c 
c 
c 


c 

c 

c 

c 


c 

c 

c 

c 

c 


area . .. surface area of semi hemi sphei e T rom theta * to bO 

degrees < m ‘ 2 ) 
rnose . - . nose radius (m> 
h wa 11 . . . wall en t h a 1 p y ( Eng 1 i sh u.n its) 

(All temperatures are in degrees Kelvin) 
tbnc-2 . temperature of fluid entering manifold from nozzle 
twmax ... maximum allowable wall temperature 

tb in ... fluid bulk temperature exiting manifold (enter ing tubes 
tbinl ... guess for tbin (used in iteration/ 1 

tfilmn, tfilm ... film temperature at which properties are evaluated 
ttube ... tube wall temperature 

no z d . . . nozzle diameter Cm) _ . , 

reyn ... Reynolds number based on nozzle diameter and nozzle exit 

velocity 

zn ... no z z 1 e— t o-t ar get separation distance (m) 
e m i ss ... r a d i a t i v e e m i s s i v i t y o f wall 
dtarg . „ . hemisphere diameter ( m ) 

qstaq, qrad , qnose ... stagnation, radiative, and nose heat 
f 1 u >: e s ( E n q 1 i s !' \ u n its) 

nuavg ... average Nusselt number based on correlation 

k ... thermal conductivity of fluid <W*m/h> 

htc ... heal transfer coefficient; of fluid < ( ' 2 > /'k > 


d l f f 
mf r no 


used in 


i t e r a t i o n 


m a s s f 1 o w r a t e a t n o z z 1 e < k g / s > 

TBN02 , TWMAX , NOZD , RNOSE , MFRNQZ , 


I r»p u t s r eq u i r ©d : 

ZN - Eli 133 

Ou t p u t s TR 1 N 

5 u t ■ , r c , n I i , ; e N A • i 1 ( c h o x c e * t b n o z ■» t w me x , n o z d ? r n o s e - m T r n oz ? z 1 1 , e m i - s , 

1 1 u b e , t b i n , e p s 1 o ri ) 

Real t brio 2 • tw.max , twmsx 1 , nozd , rnose , rnose 1 , reyn, zn,enus&, ttube , 
t b i -i p i , s i q fTi as t h e t SjS r © a s 1 1 w all . q s t a g ■» q » d * q i f o s © •> m i i i io -- v 
nuavq , t:b i. r.l , tf i 1 mn , k , h t;c , dif f * tf i lm,epslor. 

Integer cho ir e 


V a r i a b 1 e s i n d e p e i » d e n t o f f 1 u i d tern p e r a t u r e a r e cal c u 1 a t e d 
a v i d c o 1 1 v e r ted i n t o S I u \ i i t s 


a i ~ 3.1 h 1 5 9265 A 
s i q m a - 5 . 6 7 E ■- O 3 
theta ~ 25.0 * (pi/lB0.0) 


ORIGINAL PAGE IS 

OF POOR QUALfTY 


H . . . 1 11 quantities represent English units in use; note that hwall 
and qstag are both in English units. They are converted into 
m k s un i t s i n q no se b y .2931 . 

r no s e 1 = r n c« s e * ( 3 .221) 
d t a r q = r n c* se *• 2 . 0 

area "= 2.0 * pi * ( rnote**S ) * <1.0 - COS (theta) ) 
t-wmaxl - twmax * 1.8 

hwall = 0.239 * twmaxl + i.OE-05 * (twmaxl**2> 

qstag = ( ( 1 28 . 0 / < SQRT ( r nose 1 > ) > * <1-0 - < hwa 1 1 / 2080 ) ) ) 

qstaq = qstag * 3600. 

qrad = emiss * sigma * < t wmax*-*9 ) * area 

qnose = qstag ■* rnosel * 2.9816 * . Efdl. 

C o r r e 1 a t i o *n i s o n 1 y g o o d f o r u N / NO^D . 1 e . 7 . 1 J 

tfilm - 0 . 5 * ( 1 1 u b e + t b no z ) 

call d v nv < c ho ice, tf i 1 m - v isc ) 

reyn = (9.0 mfrnoz)/(pi * vise * nozd) 

i f ( < z n / no z d ) . 1 e . 7 . 0 ) t h e n 

nuavg == 2 .98 * < r e yn**0 . 585 ) * ( < no zd / d t ar g ) ** 1 . 1 0 ) * 

( ( zn/ no zd ) v —0 . 007 > > 


toy 


e I se 

wi " 1 000 ) 

1 000 format!//, ’ **** NOTE: ZN/NOZD is greater than 7****’ , 

. /,' (ZN = no z z I e— to — target separation distance? NOZD = nozzle 

.diameter ’ Please re-input data: ZN, NOZD = •’ ) 

read* ? zn , no zd 

f go to 20 
end i f 

C Now we iterate for TBIN ( temperature of the fluid exiting the 

C manifold) 

C 

tbin = 1 . Utbnoz 
10 tbinl = tbin 

t f i 1 m n — O • 5 (0.5 * (tbinl + t b n o z ) *+■ 0.5 t w rn a ; i ) 

c a 1 .1 THCOND ( c hoice? t f i 1 mn «. k > 
htc = (nuavq * k) /nosd 

t b i n = t b no z + ( ( q no se - q r ad > / < a rea * htc)) 
d i f T (A i B b ( t ta i v i — t b i n 1 ) 

C 

C Ca n c h a n q e t h e n e x t 3 i n e f o r a m ore p recise & t i m e - c o n s u \ r \ i n g 
C iteration 

c: 

i f ( cl i f f . g t . e p a 1 o n ) go to 1 0 
r e t u r n 
end 

c 

c 

c C 0 0 L A N T F 1 R 0 F E R T I E S: 

c S U B F; 0 LJ T I N E S / C U ROE FIT: 




°*iqinal p,.. 


1 


s n b r o u t i n e t h c o n d ( c h c* ice * t m a v JO 
r ea 1 k , 1 k; , a 1 < 1 0 ) * a£ < E > , a 3 ( E 1 , t ma v 
i n teg e r c h o i c (a 
d a h a a E 

* / . A h F.i r_? 4 n 3 . 1 v 29 e — 4 / 
d a t a a 3 

* / 0 . I 638 :f < -8 . 329e-4 / 

i f ( c h o ice. eq . 1 ) t h en 
a — 1 . 6030e+02 
b = -0.64379 
c = 4 . 1 l.OEe-03 
d - —7 . 679e— 06 
e = 6 . 27 79 e —09 
f = -1 . 8 19e- 12 

k: = a + (b*tmav ) + ( c*( tmav**2 ) ) + (d*( tmav**3 ) ) + (e* ( tmav**4 ) ) + 
( f * ( tmav**5 > > 
k = k * le-03 
el seif (choice.eq.2) then 
call fit4(tmavi a 2 , k ) 
e 1 se 

call f i 1 4 < t ma v i. a3 , k ) 

end i f 
return 
end 


su b r o u t i ne sp ec h t ( c: h o ice»tmavpcp) 

- r~. . 1 1 < 1 ^ I V i 1 Til '•< CO , CO \ o o i 4 - - 


lOl 


1 1 ( i. t.-‘ g t.T’ ! c r : o i c e 
data 

* /3.651 1 E03 i 2 - 536 1 E— 09 / 

if (choice.eq. 1 ) then 
a = 5.9055e+03 
bb = 1.9569E+02 
c = — 0 - 78868 
d = 1.8103e-03 
e = — 1 . 835e— 06 
f = 6.7750e-10 

Cp = a+ ( bb*tmav ) + (c* ( tmav**2 ) > + < d* < tmav**3 ) ) + ( e* ( tmav**4 ) ) + 
<f*< tmav**5> > 
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el seif (choice.eq .2) then 
i f ( tmav .It. 2 80 . ) then 
cp = 3930. 
r e t u r n 
else 

call f i 1 9 ( t ma v < a2 * c p ) 
end i f 
e 1 se 

CF' ~ 36 . 2358 ■#* ( Tmav * 0 . to 68 7 ) 

end i f 
r e t u r ri 
end 


s la b r o u t i n e d y \ t v ( c h c> ice , t m a v , m u ) 
r ea 1 mu , v i sc , a 1 ( 1 0 > - aE ( E > * a 3 ( 3 > 
i nteger cho ice 
data a B 

* / 1 <4 . 995 , -0 „ 02999 / 
data, a 3 

* / 1 . 0 1 2 e 5 > -e . 07 e— 2 / 
if ( c h o i c e . e q . 1 > t h e n 

a = 2 . 0 0 5 3 e + 0 E 
b = - 1 . 30292 
c - 3.5588e-03 
d •- —3.2 0 5 e — 0 6 
e = 9.3978e-10 

mu — a + *•. b t m a v ) + ( ( tnisv^^S ) ) + ( d ** ( t m a v * 3 > ) *+* ( e *■ ( t m a v * * 9 ) ) 
mu = mu * 1.0E-07 
r e t u r n 

elserf ( c ho ice. eq . 2 > t hen 
call f 1 1 9 ( tmav * a2 , mu ) 
return 

el se 

if ( tmav . le . 300 ) then 

mu a 10 * 2.7355e37 * (tmav ** (-16.5512)) 

else 

call f i 1 9 < tmav , a3 , mu ) 
mu = 2.679e9 * (tmav ** (-9.7321)) 
end i f 
end i f 
return 
end 


subrout i ne dense ( cho ice, tmav , ro ) 
real a 2 ( 2 ) , a 3 ( 3 ) , r o , t m a v 
integer choice 


tOf) 



¥■ 


d a t a d 

/ 1 . 171 9E3, 


142E-04/ 


if (choice.eq. 1 ) then 
a = 1 .03740+02 
b = -0.9265 
c = 3-71 98e-03 
d = -7 . 1 B9e— 06 
e = 6 - 53 1 3e-09 
f = -2.2350-12 
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ro = a+ (b*tmav)+(c*< tmav*-*2 ) ) + ( d* ( trriav*-*3 ) ) + ( e*- ( tmav**4 ) ) + 
( f * ( tmav-fr*5 ) ) 


ret urn 

elserf ( c. h o ice- eq . 2 ) then 

call f i 1 4 ( tmav?a2j r o ) 

e 1 50 

R0= 1173.0803 - (0.5093 * Tmav) - (2.795E-04 * (Tmav ** 2)) 


end i f 
r e t u r n 
end 


5 i t b r o u t i n e f i 1 4 < t m a v ? £< ? a n s w e r ) 
rea 1 tmav ? a ( 2 ) •» answer 

answer = a ( 1 ) * (2.71 82 S 1 828 ** ( a ( 2 ) * tmav ) ) 

r e tur n 
end 


/o<? 
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Coc lant considered : Liqu id hydrogen 
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Appendix 3.2: Ethylene Glycol Design 
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Appendix 4: Hydrogen properties 
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Table 4.12 (continued) 
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Appendix 5: Hypersonic Cruise Vehicle Configuration 


ijypersonic Cruise Vehicle fHCV) Confinuratinn ( \Mrijten \>y Lvra 
It was determined early in Winter Quarter 1988 that the 
Advanced Aeronautics Design Group (AAD) needed a unifying concept for the 
design of the hypersonic drone. This was accomplished by choosing a 
configuration that had been researched extensively at NASA Langley 
Research Center in the late 1970’s. The advantage of a single tested 
configuration for the entire class is clear: each group would base their 

designs on the same vehicle, with the same restrictions applying to all 
three groups, thereby forcing greater interaction among the groups. 

The configuration chosen by the AAD students was a 70 degree 
swept delta wing design conceived by NASA Langley engineers. It was a 
hypersonic research airplane concept, and has been dubbed the Hypersonic 
Cruise Vehicle (HCV) by the AAD students. Figure 1 shows the HCV as it 
was originally conceived. The HCV was designed to be air-launched from a 
B-52 at Mach 0.8 and rocket-accelerated up to a cruise speed of Mach 7 by 
four liquid hydrogen rocket motors mounted in the base. At this point, the 
scramjets^ would be ignited and the vehicle would cruise until burn-out, at 
which point it would return to base for an unpowered landing. A detailed 
description of the HCV’s configuration is given in Figure 2. In developing 
the HCV design, a wind tunnel model was built, and performance 
characteristics of the model were collected and tabulated in Reference 1. 

The geometric characteristics of the model are tabulated in Table 1. The 
scaling for the model was 0.021 of the size conceived for the actual HCV. 

The first step in adapting the HCV for the purpose of a 
hypersonic drone was to size it according to the needs of the class. This 
included calculating the amount of fuel required to go from the air-launch 
speed to cruise speed. From that, size and weight necessary for the fuel 
requirements were determined . In order to do this, it became necessary to 
determine a means of air-launch. The B-52 has size and weight 
constraints consisting of a 24-ft. wing span limit, a 9-foot fuselage 



height limit, and a gross weight restriction of 70,000 pounds. This limits 
the size of the drone considerably and, in order to reach Mach 10, it was 
found that the weight limit was exceeded. Alternatives were then 
examined. The 747 and B-1 aircraft were determined to be viable options 
for the air-launch. 

Therefore, the HCV is a basis upon which each group can build 
their design. It gives the Flight Systems group the geometry and spacing 
of the vehicle so that they can determine which sensors and how many are 
required in certain places. This would result in a more detailed and more 
exact design. The Thermal Management group can determine critical areas 
on the HCV that would require cooling. They can change the design within 
certain reasonable limits to satisfy their cooling system design. A* 
definite configuration for the the drone gives the Propulsion group a 
starting point on the design of the underbelly for the vehicle; that is, a 
first iteration could be performed from which they refined their design to 
achieve an optimum configuration. 
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Introduction 


Composing on© third of the UCLA AAD program is the Flight 
Systems Group. The purpose of this group is to design, find, and 
adapt the necessary sensors and determine methods of integrating 
these sensors with the structure of the proposed drone envisioned by 
the class. The sensors to be designed are to complement the unique 
needs of the Propulsion and Thermal Management Groups. 

The hypersonic flight regime in which the drone is designed to 
be flown and tested makes the design of the vehicle and its systems 
particularly challenging. The extremely high speed environment will 
lead to a host of problems, including high temperatures, pressures, 
and static and dynamic loads. From the viewpoint of sensor 
systems, many of the classic methods are inadequate in this region. 
For example, because of the huge heat fluxes, external probes are not 
usually feasible. Instead, non-intrusive methods need to be 
examined. 

The parameters to be measured that were studied by the group 
were determined by the needs of the other two groups. These 
include temperature, pressure, flow constituents and 
concentrations, velocity, and skin friction. Sensors were designed 
from such diverse fields as quantum electronics, semiconductors, 
and fiber optics. The sensors researched originated from a variety 
of environments, from medical laboratories, to experimental 
laboratories, and our own class members. 

This report will now detail the various sensors researched 
developed, and adapted by the Flight Systems Group for the 
hypersonic drone design. The report begins with the very flexible 
and powerful sensor system called CARS. Next are two conceptual 
designs for shock-wave and pressure sensors, followed by a fiber 
optic pressure sensor. A hydrogen-concentration sensor realized 
using a modified semiconductor device, and a device using a sapphire 
diaphragm. The final sensor detailed is a laboratory sensor for 

measuring skin friction such as that encountered in the combustor of 
the scramjet. 
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Coherent Anti-Stokes Raman Spectroscopy 


Coherent Anti-Stokes spectroscopy, or CARS, is a laser 
technique which is capable of measuring velocity, temperature, and 
constituents of gaseous flow. This is done by measuring the 
returning characteristics of emitted laser light from the molecules. 

The range of these parameters (velocity, temperature, 
pressure, and species concentration) vary under specific conditions. 
Available data reveal as follows: 

A. Temperature and Pressure (in N2 gas) 

In N2 gas, measurement error is estimated to be less 
than 1.5% in time-averaged measurements, and from 3% to 6% for 
instantaneous measurements, over a wide range of temperature and 
pressure. 

B. Species Concentration 

Species concentration sensitivity is greater than 1% to 
5%, over a limited range of concentration (0.5% to 20.% depending 
upon the species). Usually, this is sufficient for practical 
combustion systems. 

C. Velocity 

Velocity of a flow constituent, or species, is determined 
from the Doppler shift of the returning beam, and hence, is subject 
to the limitations in range of detection possible for this type of 
frequency shift, if any. 

D. Pressure 

Though pressure is probably better dealt with using 
Inverse Raman Spectroscopy (IRS), a CARS system may be capable of 
measuring variations of pressure in a flow. 

The accuracy of the CARS measurements depends largely on the 
prevailing conditions, the type of CARS used, and the method of 
processing the CARS data. However, good accuracy of readings under 
conditions of interest to the Flight Systems Group have been 
obtained: 
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"The velocity in a supersonic jet flow has been measured using 
Continuous Wave CARS and agrees to within 4 percent of the 
calculated value." (ref. 8). 

The dimensions of the CARS apparatus, i.e. its size and weight, 
depend upon the type of laser used, as well as the geometry of the 
optics technique used to generate and process the CARS signal. The 
laser can be made fairly small, thus, it would not have as much 
influence on the size of the CARS set-up as the rest of the 
apparatus. A reasonable estimate for the area of the entire set-up 3 
ft by 2.5 ft. The entire system size does not seem to be a rigid 
factor, however, since laser technology is moving toward smaller, 
and smaller, optical equipment. 

The weight should not be more than a few pounds for the laser, 
and a little more for the optics that direct the beam, excluding any 
anti-vibrational equipment necessary for the implementation of any 
sort of laser measurement system on board the drone. 

Power consumption (excluding the computer analysis power 
drain) would depend upon the number and type of lasers used. The 
supersonic jet flow measurements mentioned previously were taken 
with a 100 mW laser. 

The CARS measurement rate varies with the type of CARS, and 
the laser pulse frequency used. At least one type of CARS can be 
used for continuous measurement. This requires a continuous-wave 
laser instead of a pulsed laser. 

Whatever type of laser is used would definitely require an 
active cooling system, as very high power would be used. 

Finally, we consider the placement and distribution of the 
system. The laser would be inside of the drone in a cooled area. For 
measuring wide areas of flow, a rotating platform would enable one 
CARS set-up to suffice. Or, alternatively, more than one clear 
window could be installed, and the emerging CARS beam re-directed 
through any of these temperature-resistant quartz windows. (Fig. 1 ) 

CARS has several advantages that make it a good sensor for 
the drone. For instance, CARS will perform a few types of 
measurements at once, is generally non-intrusive, and the beam is 
not damaged by intensive heat. It is possible that CARS need not 
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have any apparatus to reflect the beam back from the medium it 
measures. This is due to the fact that the molecules in the medium, 
themselves, generate a laser-like beam. 

The limitations of the CARS technique depends on which type 
of CARS is used. Some types work well in high pressure, turbulent 
flow, while others will not work under those same conditions. Most 
of the major limitations of CARS are a result of pressure gradient 
problems (the beam is refracted by pressure gradients, resulting in 
"beam steering" problems). 

CARS will satisfy the needs and specifications of the 
Propulsion and Thermal Management groups. For example, in the air 
flow boundary layer outside of the drone, refraction and reflection 
("beam steering") could possibly be used to our advantage. This is 
because the return beam is reflected/refracted back exactly in the 
same path. CARS spectra, and some geometry, can give us 
temperature and boundary layer data. 

CARS can also measure flow constituents in the combustor 
nozzle, the mixing rate and turbulence of the fuel, and the 
concentration of elements downstream and upstream of ignition. 

The most important aspect of CARS is the Raman effect. Laser 
light has a single frequency and phase quality. When light of a 
particular frequency hits a molecule, it emits a laser-like beam 
with the same qualities of the original beam, except that the 
emitted light is shifted by a frequency that is unique to that 
particular molecule or atom. This phenomenon is known as the 
Raman effect. 

When this light hits another molecule of the same type, the 
probability of the molecule emitting a Raman shifted beam, 
originally low, greatly increases. The effect is amplified as the 
beam traverses through a sample of these molecules. 

The beam is returned either by a reflecting apparatus at the 
other end, or by some technique that directs the beam back towards 
the source. The intensity, phase, and rotational direction of the 
beam and its spectrum gives information regarding temperature, 
concentration, and pressure of the molecules hit by the laser. Most 
of this information is determined by comparing the beam spectrum 
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to some theoretical or reference spectrum. Velocity is measured by 
noting the Doppler shift in the returning beam frequency. 

Coherent Anti-Stokes Raman Spectroscopy, or CARS, is a laser 
spectroscopy technique. Spectroscopy is the study of spectra 
obtained from a spectroscope. A spectrocope is an apparatus that 
takes incoming light and by using a prism or diffraction grating 
forms a spectrum by separating the various wavelengths of the light 
into bands of varying intensity. The wavelength of a particular band 
can be determined from calculations using the position where the 
band occurs in the spectrum. How well these bands are 
distinguished from one another determines the resolution of the 
spectrum. 

In CARS, two or more lasers are combined using a method 
called Optical Wave-Mixing to generate a single laser beam. This 
beam has a frequency that molecules in a media will resonate when 
hit. The enhanced beam that results from the molecule is the CARS 
beam or signal. This resonant response is the result of the Raman 
effect discussed earlier. Spectral analysis of the returning CARS 
signal can reveal information about the temperature, species 
concentration, pressure, and velocity of the molecules that it 
originates from. 

The two or more beams that are optically wave-mixed are 
subject to certain requirements. For example, in two-laser CARS 
one laser is chosen to be the "pump" laser. The other laser is known 
as the "Stokes" laser because its frequency is shifted in the "red" 
direction of the electromagnetic spectrum with respect to the pump 
laser. The Stokes frequency is not fixed, and can take on different 
values. 

The wave-mixed beam goes through the medium to be examined 
and results in a CARS beam and frequency. This is the "Anti-Stokes" 
beam, because its frequency is twice the pump frequency less the 
Stokes frequency. It is shifted towards the blue end of the 
electromagnetic spectrum with respect to the pump beam. 

There are two types of Stokes lasers used in CARS methods: 
narrowband and broadband. These two ways of generating the Stokes 
frequency separates CARS into two groups which are named 
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according to the laser used. Narrowband (also known as "scanning") 
CARS builds the CARS spectra piece by piece by "tuning" the laser to 
a small range of Stokes frequencies, taking a measurement for each 
different frequency. This type of CARS is generally restricted to 
steady-state flow. 

Broadband CARS, in contrast, will generate a complete 
spectrum for a given range of the Stokes frequency in one 
measurement. While a narrowband laser will output only a very 
narrow range of frequencies, a broadband laser outputs a wide range 
of frequencies. This type of CARS is necessary when measuring 
"randomly time-varying flows", since the sample could change from 
one point in time to the next. 

Both scanning and broadband CARS usually examine only one 
species in a sample at a given time. There is a method that is 
capable of measuring all of the major constituents in a flow 
simultaneously. It is called Dual Broadband CARS. This method 
seems to hold the most promise to fulfill the sensor requirements 
for the hypersonic drone. 

Dual Broadband CARS is one of a group of CARS techniques that 
utilize more than one pump and/or Stokes lasers to create the CARS 
beam. These are multi-frequency (multi-color) methods, since each 
laser operates at its own frequency. This method also contains the 
most potential for high pressure measurements which are expected 
for the combustor of the drone. 

Alan C. Eckbreth and Torger J. Anderson, of United Technologies 
Research Center developed Dual Broadband CARS. "The technique 
employs two independent broadband Stokes sources in combination 
with a pump laser. In addition to the two-color wave-mixing 
processes between the pump and Stokes lasers, spectrally resolved 
CARS is produced in a three-color process from species whose 
Raman resonances correspond to the frequency difference of the two 
broadband sources. CARS is thus derived from a large number of 
species simultaneously removing the nominal limitations of CARS to 
interrogate only one constituent at a time." (Ref. 4) 

For the Flight Systems Group purposes, the process of two and 
three color wave mixing is not so important as is the fact that two 
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broadband Stokes sources are used with a pump laser. The result is 
the capability to measure more than one species in a sample at a 
time. 

In the article "Multi-Color CARS for Simultaneous 
Measurements of Multiple Combustion Species", Eckbreth and 
Anderson outline three other multi-color techniques as well as 
rotational/vibrational CARS methods. They are Dual Broadband, 

Dual Stokes, Dual Pump, Dual Pump-Stokes CARS, Pure Rotational 
and Combined Rotational/vibrational CARS . 

Dual Stokes and Dual Pump CARS are used primarily for 
measuring two flow constituents at once. Measuring three or more 
are theoretically possible, but the application is too complicated to 
be practical. Both of these approaches result in simpler CARS signal 
detection schemes. 

Dual Pump-Stokes CARS is interesting because the Stokes 
bandwidth can be made very narrow. In hydrogen air combustion, 
this method can then "track the appearance of H20 product in the 
reaction". This means that it could keep track of how much fuel is 
burned in the drones’ scramjets. 

Pure Rotational CARS is known for its high accuracy in 
measuring temperatures below 1000 K, even at high pressures. This 
method is difficult to apply. It uses the rotational approach on Dual 
Broadband CARS, making it simpler to implement. The advantages of 
this is the elimination of a Stokes laser, no calibration for the pump 
laser, and the pump frequency can be arbitrarily chosen for maximum 
power output. 

Combined Vibrational/Rotational CARS combines Vibrational 
CARS with Rotational CARS. The result is a method that will 
measure two constituents at once over a very wide range of 
temperatures. 

After considering the assets of all the CARS techniques, Dual 
Broadband CARS looks like the most promising as a sensor for the 
hypersonic drone. The major drawback is that there will be a loss of 
signal intensity in the CARS beam as compared to "regular" CARS. 
"Due to the increase of the CARS signal with increasing pressure, 
Dual Broadband CARS should be capable of single pulse 
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measurements of the major constituents at elevated pressures 
typical of gas turbines, internal combustion engines, and burning 
propellants. Furthermore, in Dual Broadband CARS, a beneficial 
spectral averaging may occur resulting in improved single pulse 
spectral quality". It is also noted that "dual Broadband CARS may be 
restricted, particularly for dynamic measurements, to elevated 
pressure combustion applications. Signal decrease accompanying 
broadband sources is compensated by ...scaling... with density". (Ref 
5 , 6 ) 

This suggests that the use of Dual Broadband CARS in the high 
pressure scramjet engines to monitor the temperature, velocity, and 
species concentration of the flow there. A remote sensing system is 
possible for measuring velocity, pressure, and temperature outside 
of the drone as well. 

To summarize, several CARS diagrams have been included. One 
of these is a "typical" CARS system, and also a configuration for 
measuring the after-burn of jet engine exhausts. After primarily 
examining the needs of the Propulsion Group , Dual Broadband CARS 
was selected as the method most capable of fulfilling them. 
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CARS Supporting Optics System 

Having decided upon which type of CARS to implement, the 
manner which this was to be done remained to be developed. The 
following optics arrangement was designed to make the most of the 
CARS laser spectroscopy system as applied in the scramjet 
combustor. (Fig 2.) 

The first is a result of our desire to check flow patterns in the 
scramjet in several areas. The optics arrangement is capable of 
achieving this goal because it uses the idea of one CARS system 
directed out of many windows. 

The second way the optics arrangement optimizes the CARS 
system is by making the most of the available laser power coming 
from the optically wave-mixed CARS beam. This beam is split up 
between these windows in terms of fractions of the original beam 
pulse(s), as opposed to dividing a single pulse of the beam between 
all of the windows simultaneously. 

The third way is that the CARS system is physically connected 
to a single, outer structure. This minimizes the effect of the 
aircrafts vibrations on the data , since all of the components of the 
optics system will vibrate relative to the drone and not each other. 
Further, the system may operate at such high speeds ( picosecond 
pulses, for example ) that , in comparison, any aircraft-related 
vibration would be insignificant. More precisely, if pulses are 

generated fast enough, data could possibly be gathered between one 
vibration and the next. 

Thus, we hope to achieve a system capable of mapping at least 
part of the scramjet combustor flow, with a minimum amount of 
power loss, and in the most accurate way possible. 

The way the CARS supporting optics system ( CARS S.O.S.) 
operates is rather complicated. The CARS beam is generated in a 
separate system that is housed inside the main structure of the 
CARS S.O.S. As this beam comes out of the generating system, it 
goes through an aperture (marked "A" in the figure) and into the 
Supporting Optics System. As it travels horizontally into the 
system, the beam is directed to the first rotating shaft beam 
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splitter (marked "C") by a totally reflecting mirror ("B"). On this 
shaft, and on the second , tilted axis shaft below it, are disks 
which both allow the transmission and reflection of the beam. 

These disks are designed such that the beam will reflect off each 
one , in succession, for a predetermined pulse length, L. This pulse 
length is determined by the amount of reflecting disk that is cut 
away in the disks between it and the laser aperture. For a cut away 
length of circumference, d , assuming the laser beam traces this 
circumference, a pulse will be generated of length, 

L= ( c * d )/(Pi * r *f) 

where: 

r= distance of beam from disk’s center 

c- the speed of light in a vacuum 

f= the frequency of rotation of the shaft about its axis. 

Note that the pulse length is dependent on the distance of the 
beam from the center of disk and not necessarily on the radius of 
the disk. The radius of the disk only determines the maximum r. This 
makes the optical design quite a tricky geometry problem, since 
the radius of a disk depends on the position of the disk on its shaft 
relative to the incoming beam, the horizontal spacing of the 
windows, and the amount of room between the horizontal and tilted 
axis rotating shafts. Thus, the pulse length may be changed by 
change in axis frequency, but the other parameters are fixed by 
geometrical requirements. 

Once the CARS beam is split into a pulse by the horizontal 
shaft arrangement of disks, the pulse travels to its reflecting 
horizontal disk which directs it to the appropriate second, tilted 
axis shaft disk. These tilted axis shaft disks, also, alternately 
reflect and pass the beam pulse -- only these reflect the pulse down 
through a collimating lens ( of 300 mm in focal length ) and through 
one of the temperature - resistant quartz windows into the 
scramjet combustor. There, a CARS signal is generated by the 
specific, excited molecules, and travels directly up to the tilted 
shaft disk which first reflected it down. As soon as the pulse is 
over and the CARS laser is blocked from the reflecting horizontal 
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disk s view, the tilted shaft disk allows the beam to travel straight 
up to an optical detector. 

If the pulse is long, since the CARS stimulated molecule beam 
is approximately instantaneously generated, an additional detector 
may be placed around the aperture of the original beam to collect 
any overlapping signal that occurs during the pulse generation. Any 
overlap could then be routed to the appropriate spectral analysis 
channel. It is possible to use this same system with an entirely 
reflecting second shaft and one photodetector at the CARS beam 
aperture . However, having separate photodetectors for each pulse 
has the possible advantage of easier signal separation at high pulse 
frequencies and faster electronic processing of the almost 
continuous wave of CARS returning pulses. 

Note that only one pulse will be reflected up from the 
scramjet at any moment, in either case. Hence, one spectrograph and 
computer, operating continuously, could receive all of the data. It 
may, however, take more than one computer to store .analyze, and 
process the information at the speed required. 

Note, also, that other means could have been used to make the 
division of the CARS beam between the windows much simpler, but 
with a considerable power loss increase. Beam splitters which 
partially reflect and transmit simultaneously, could have been used 
to split the original beam. But such beam splitters, used twice in 

succession, would have reduced the outgoing beam intensity to 25 % 
of the original CARS beam. 

However, in return for conserving power and the other 
optimization inherent in its design, the proposed Supporting Optics 
System has its own particular drawbacks. First, the size of the 
optics system in very much dependent on the spacing of the windows 
and number of windows used. A fixed amount of space is necessary 
to contain the spinning disks and to make sure that the two spinning 
shafts do not collide with one another. 

Also, more power is consumed than in a stationary optics set- 
up, since power must be supplied to the shaft rotary motors. 

Further, the very fast rotation of the shafts would generate heat, 
which would add to the overall problem of cooling the CARS. 
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Lastly, such a system would be more prone to mechanical 
failure than a stationary system using half-reflecting/ half- 
transmitting mirrors. 

In summary, however, even with these drawbacks (increased 
power consumption, greater cooling requirements, and possible 
mechanical failure), the proposed optics arrangement for the 
application of CARS in the drone's scramjet combustor should be 
able fulfill the requirements it was designed to to a satisfactory 
degree. It will allow the measurement of combustor flow 
characteristics in several different areas in the combustor, provide 
the maximum possible power to these measurements ( and, hence, 
the more accurate results ), and minimize the effects of any 
aircraft-related vibration. 
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The next two devices described are different than the other 
sensors. Both use novel pulse-laser techniques that were designed 
conceptually from scratch. Thus, there is considerable ambiguity as 
to the actual feasibility of the devices since neither the dollar cost 
nor any available substitutes have been considered. Design is 
strictly conceptual that can be used as a springboard for further 
research. 


An Optical Pressure Sensor 

This device (shown in Fig. 3) uses a Michelson-type 
interferometric technique to measure the displacement of a 
material, such as quartz or graphite, under the load of substantial 
pressure. The sensor does so by an arrangement of thin aluminum 
mirrors of proper reflectivities embedded in the bulk material 
through which a laser beam is launched. Fringes are formed as 
pressure displaces the upper mirror downward towards the laser 
source. These fringes are then counted and the measure of bulk 
material displacement under the external load is then calculated. 
Through accurate calibration or theory, the displacement of the 
mirror, which moves in parallel with the bulk material it is 
embedded in, is related to the pressure the bulk is exposed to. 
Therefore such a device is useful in measuring pressures exceeding 
ordinary pressure sensors, pressure as high as the bulk material can 
withstand while still retaining its elastic properties. Outside 
dynamic pressure anywhere on the surface of a hypersonic airplane 
can be measured with such a device. The data can be sampled at 
rates above 100 KHz (limited only by the propagation delays of the 
associated electronics). 

This device suffers from noise that results from all mirrors 
being displaced simultaneously from pressure. One way of 
correcting this noise-like effect is by inserting a soft buffer layer 
between the upper mirror and the lower part of the device. 
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Multiple Shock-Wave Sensor 


The operation of the device is dependent of the shock wave 
creating a boundary layer. The differing physical properties of 
abrupt pressure and temperature changes across the boundary layer 
leading to a spatial discontinuity in the index of refraction. A laser 
beam can be reflected off of the shock wave because of the index 
discontinuity. 

Due to the high hypersonic speeds, the shock waves will "hug" 
the plane, being nearly horizontal with respect to the surface. A 
pulsed laser beam is launched upward at right angles to the surface 
of the drone into the air (as shown in Fig. 4). The laser beam will be 
partially reflected upon each impact with each shock wave that it 
encounters. By detecting the back-scattered beam, its intensity- 
time profile can be used to predict the locations of the shock waves 
away from the plane. Knowledge of the appropriate permeability 
coefficient is needed in the calculation of the shock wave distance. 
This can be easily obtained with a simple in-flight measurement. 

Noise from various sources such as the beam oscillating 
between two shock waves and beam-steering may be mis- 
interpreted for a shock wave. The high temperatures and the near- 
plasma environment outside make the device difficult to design. It 
is difficult to single out only the intensity peaks produced by 
scattering the beam from a shock wave front, effectively rejecting 
noise contributions. Beam steering effect can lead to significant 
errors in the measurements due to the turbulent environment 
expected. 

The velocity of the drone can be conceivably predicted from 
the angle the shock wave makes with the fuselage. By scanning the 
shock wave with the arrangement described above, such an angle can 
be computed. Combined with the altitude of the drone, and other 
basic parameters, the velocity can be predicted. 
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Fiber Optic Pressure Sensor 


Fiber optics are currently being used in the measurement of 
many quantities, such as temperature, acceleration, velocity, 
electric and/or magnetic fields, fluid levels and acoustic waves. In 
the measurement of these parameters, the fiber plays either an 
active or passive role; active when the fiber itself is the sensing 
element and passive when the fiber simply transmits information 
from the sensing element. The device considered here is a pressure 
sensor in which the fiber takes an active role in the measurement. 

The sensor is a flattened hollow cylinder about which is 
wrapped an optical fiber. At one end of the fiber is a source, an 
AIGaAs LED with a power consumption of the order of milliwatts. At 
the other end is a detector, a Si PIN photodiode, the pressure 
differential between the inside of the cylindrical cavity and the 
outside of it causes the radius of curvature of the flattened section 
to decrease from infinity to some finite value. As this radius 
decreases further, there is a corresponding approximately 
exponential increases in the attenuation of the intensity of light 
that reaches the detector. The light attenuation at the curved ends 
is constant and is taken into account when calibrating the device. 

The amount of light reaching the detector is therefore a function of 
the pressure difference between the inside and outside of the cavity. 
The outside of the cavity could be another cavity at a reference 
pressure, or a cavity in which the pressure can be varied to increase 
the range of measurement of the device. See figure 5 for a picture 
of the fiber optic pressure sensor, (ref. 12) 

This sensor can be adapted to a high temperature environment 
by constructing the cylindrical cavity from titanium and threading 
the approximately 100 micron diameter fibers through tubes welded 
to the outside of the cylinder. Water or the liquid hydrogen on board 
would be pumped through these tubes to keep the fibers at a 
relatively constant temperature. In this way the device can be 
located near the surface of the inside walls of the engine or the nose 
of the craft for a flush air data system to optimize the pressure 
propagation delay time. The source and detector would be located 
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away from the high temperatures, fibers of course being used to 
transmit light to the sensor. In this way, real time measurement of 
pressures on the surface of the aircraft or inside the engine could be 
made. 

The theoretical sensitivity of this device is given by the 
following equation: 

S*(1/l)(dl/dP)= Sopt (1- ji^/ehS) r 2 (fl + i2/7i)/(IB + 7 t) 
where l=intensity reaching detector 
p=pressure 

Sopt ■ normalized optical sensitivity for a 
constant length L of fiber bent to radius R. (ref. 7) 

This device will use single-mode fibers rather than multi- 
mode fibers which were used in the Wright experiment (ref. 20) 
because of the vibrations of the multi-mode fibers encountered 
during flight would effect the electromagnetic power propagation 
and thus pressure measurement. Two advantages are gained, 
however, with the use of single-mode fibers. First the sensitivity, 
S 0 pt is an order of magnitude greater, 20 compared to 2. Secondly, 
with single mode-fibers is proportional to the number of half 
turns the fiber makes around the flattened cylinder, so that the use 
of titanium, with two orders of magnitude greater Young's modulus 
than the epoxy resin of the Wright experiment, causing a two orders 
of magnitude decrease in the sensitivity can be compensated for 
with increased half turns of the fiber. Since is already an order of 
magnitude larger, two extra half turns brings the sensitivity of this 
device into the range of the Wright experiment, S=8*/bar (1 bar = 

100 kPa). 

The dynamic range of the epoxy resin cylinder was tested for 
200 kPa. This range could, of course, be extended with the use of 
the titanium cylinder. 
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MOSFET Hydrogen Concentration Senducer 


Hydrogen flow concentration inside the combustor of the 
SCRAMJET is a very important parameter to be measured because it 
is related to the rate of combustion, the energy consumption, and 
the engine efficiency. The sensors used in the are relatively big and 
intrusive. The MOSFET hydrogen concentration senducer (ref. 3) has 
many advantages such as small size, low power consumption, high 
sensitivity, fast response, and relatively non-intrusiveness. 

A senducer is a sensor and a transducer combination which 
senses physical quantities and then converts them through chemical 
parameters into electrical signals. The acronym, MOSFET, stands for 
Metal-Oxide-Semiconductor Field-Effect Transistor. It is a four- 
terminal device as shown in Figure 6. 

The operation of a MOSFET is controlled by a vertical action and a 
horizontal action. The gate-bulk potential difference Vqb sets up a 
capacitive effect which adjusts the amount of carriers (in our case, 
electrons) in the conducting channel. There will be carriers in the 
channel if, and only if, the gate-bulk potential is greater than, or 
equal to, the threshold voltage Vy. This is the vertical action. The 
horizontal action is controlled by the potential difference applied 
between the drain and the source VdS- Current flows from the drain 
to the source if vertical action causes carriers to be present in the 
conducting channel. Depending on the value of Vds. there are two 
regions of operation--the linear region and the saturated region. 

The linear region occurs if Vds £. Vgs - Vj and the saturated region 
occurs if Vds > VGS - Vj. The drain current equations are given by: 
ID - ji( W/L)C0X [VDS (Vgs - VJ) - VdS 2 /2] linear region 
ID = (1/2)(W/L)Cox (VGS * Vj)2 saturated region 

Where p. » mobility of electrons 

W/L - gate width / gate length 
COX - permeability / thickness of Si02 
Notice that the current in the linear region is less than the 
current in the saturated region and thus a MOSFET operating in the 
linear region consumes less power. Also, the fact that the current 
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in the linear region depends on Vps makes the linear region more 
flexible for out design. 

In our design, the gate material used is palladium (Pd) instead of 
polysilicon used in microelectronic devices. Palladium, when in 
contact with hydrogen flow, acts as a catalyst for H2 ionization t 
give 2H+. This H + diffuses into the thin Si02 layer to give a charge 
difference AQ. A AQ changes the threshold voltage Vj. If at AQ = 0, 
= ^TOi then Vj - Vjo - Vj ■ AQ/COX- Now, if we substitute Vj = 
Vjo - AQ/Cox into the linear region Id equation, we have: 

ID = P(W/L)C0X[VDSAQ/C0X + VDS(VGS - VJO - Vds/2)] 

Notice that for a particular Vds, Id is linearly related to AQ with a 
slope equal to |i(W/L)Vds and the plot is shown in Figure 7. 

The senducers can be embedded in the walls of the combustor. 

The advantage is that although the sensor cannot stand high 
temperatures, we need not worry about heating effects because the 
combustor walls will be actively cooled regardless of the presence 
of the MOSFET senducers. 


Note that we can use platinum instead of palladium as the gate 
material. However, since they act the same for our purpose and that 
platinum is much more expensive than palladium, palladium is 
obviously a better choice. 


The senducer has the advantage of small size— approximately 
10 -15 m 3 f | 0W p 0wer consumption on the order of 0.1 mW, fast 
response time on the order of a few nanoseconds, and high 
sensitivity to H2 concentration variation because molecular 
relaxation takes only a fraction of a nanosecond. 


To ensure accuracy, each senducer should be calibrated 
individually because of the doping variations and difference in 
threshold voltages. However, if all the senducers are made on th,e 
same silicon wafer, this calibration problem can be greatly 
simplified. 
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Sapphire Diaphragm for Pressure 
and Temperature Measurement 

A sapphire diaphragm approximately an inch in diameter is 
placed flush with the aircraft's skin. Pressure measurement with 
this device is based on the photoelastic effect. Light from an LED 
source is linearly polarized 45 degrees relative to the surface of the 
diaphragm and passed across the diameter of it. As pressure is 
applied to the surface, the linearly polarized light becomes 
elliptically polarized due to the photo-elastic effect. A second 
polarizer , again placed at 45 degrees relative to the surface, filters 
out the elliptically polarized light at the opposite end of the 
diaphragm. The degree to which the linearly polarized light is 
attenuated is thus a function of pressure. The accuracy of this 
device can be increased by using two LED sources and alternately 
polarizing them +45 degrees and -45 degrees relative to the surface, 
therefore, obtaining four outputs. The resulting four outputs can 
then be averaged. Temperature measurement with this device is 
based on blackbody radiation. As the temperature of a substance 
increases the amount of radiation emitted increases and the 
wavelength at which it is emitted decreases. The relation that 
describes this phenomena is Planck's equation, (ref. 17,18) 

Since sapphire is not an ideal blackbody, the above 
relation must be multiplied by the emissivity of sapphire, 0.15, 
which is an indication of the degree to which it approaches the ideal 
blackbody. 

Sapphire was the material of choice for several reasons. First 
and foremost is its high temperature properties. Its melting point 
is 2030 degrees Celsius with its maximum normal use temperature 
at 1950 degrees Celsius. It is unaffected by oxidizing and reducing 
atmospheres at elevated temperatures. Its crystalline structure, 
trigonal, is identical to that of alpha-quartz which was used in the 
photoelastic experiments. The tested range of quartz was 10 MPa 
with a 0.07% accuracy. There is no reason to believe that sapphire 
would be any less. It is a low power consumption device; the LED's 
are lOOmW sources. And finally, since the diaphragm is flush with 


19 - 



the surface there is no need to wait for the pressure wave to 
propagate down a tube to the pressure sensor located inside the 
aircraft. 
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The Dual-Laser-Beam Skin 
Friction Interferometer 

The skin friction sensor, using a non-intrusive double-laser- 
beam oil viscosity balance technique, was developed by Monson and 
Higuchi (ref. 14). Two out-of-phase laser beams, reflected from an 
oil film, gives alternative constructive and destructive 
interferences. The number of fringes, and thus the rate of change of 
the oil film over a given time span, can be precisely computed in 
terms of the known laser wavelength. The skin friction is related to 
the rate of change of thickness of the flowing oil film. The oil film 
used is silicon oil since it has a wide range of viscosity, a relative 
insensitiveness to temperature, a low surface tension, and a very 
low vapor pressure. The following paragraphs give a brief outline of 
its operation and characteristics. 

Figure 8 shows the configuration of this device. A He-Ne laser 
beam is generated. The laser beam is expanded by a telescope, 
linearly polarized, and passed through a neutral density filter to 
reduce the beam power level to a value that avoids excessive oil 
heating. The beam is then split into two, using a interferometric 
flat. One of the beams is rotated by 90° by the half-wave 
retardation plate, and both beams are focused on the oil film. The 
reflected beams are separated by a polarization beam splitter and 
focused on a photodiode after passing them through polarizers and 
narrow band interference filters at the detectors for additional 
noise reductions. The interference patterns are then analyzed to 
give skin friction. 

This device is very accurate and inexpensive. It can be used in 
three-dimensional laminar or turbulent flow. It is only slightly 
intrusive and can be used at various points on the combustor surface. 

The only problem occurs when there are dust particles in the oil. 

A dust particle exactly at a beam focal point could cause erratic 
interferometer behavior. A large dust particle just ahead of a beam 
focal point could locally perturb either the oil, or the oil flow, 
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enough to influence the skin friction measurement, but could not be 
detected on the interferometer record. 
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Conclusion 

This ends the report for the Flight Systems Group. The sensor 
system design was approached with an integration with the 
Propulsion and Thermal Management groups. The challenging 
environment and the variety of needed measurements called for a 
equally large variety of sensors from differing backgrounds and 
applications. These innovative designs should provide an interesting 
testbed for the sensor systems on the hypersonic drone. 
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